Esta proyecto corresponde a la PRAC_1 de la asignatura Minera de Datos del Master de Data Science de la UOC (Universidad Oberta de Cataluña)
La idea inicial era crear mi propio dataset pero he cambiado de opinión debido a que iba a dedicar más tiempo a recopilar todos los datos y crear el dataset que el propio modelo en sí.
Por ello he decidido investigar las websites propuestas para la PRAC y de entre todos los dataset disponibles me he decidido por el existente en UCI Machine Learning Repositorty , “Heart Disease Data Set”.
El termino cardiopatía se refiere a cualquier enfermedad del corazón independientemente de su naturaleza . Según la OMS las Enfermedades cardiovasculares (ECV) son un conjunto de trastornos al corazón y de los vasos sanguíneos, siendo la principal causa de defunción en el mundo.
Como datos llamativos:
• 17,5 millones de personas murieron por enfermedades cardiovasculares en el mundo.
• 80% de los infartos de miocardio y de los AVC prematuros son prevenibles.
• Más del 75% de las muertes causadas por ECV se producen en países de ingresos medios y bajos.
En España y como podemos leer en la web de la Sociedad Española de Cardiología (SEC): La enfermedad cardiovascular continúa siendo la primera causa de muerte en nuestro país, con 120.859 fallecimientos registrados al año, según los últimos datos del Instituto Nacional de Estadística (2018).
Si nos fijamos en el segundo punto de las afirmaciones de la OMS, el 80% de los infartos se pueden prevenir. Esto significa que con un estudio pormenorizado de los factores de riesgo, utilizando la cantidad de datos médicos existentes y modelos de machine learning, estas causas de defunción se podrían disminuir considerablemente pudiendo prever con antelación que grupos de población con sus respectivas características médicas van a tener más riesgos que otros.
Por ello, podemos considerar la predicción de enfermedades cardiovasculares como uno de los objetivos más importantes de la mineria de datos y del machine learning dentro de la medicina.
Por otro lado, este dataset y su estructura nos permitirá realizar modelos de clasificación y clustering, modelos aprendidos en la asignatura de minería de datos cuyo proyecto nos permitirá afianzar lo aprendido.
Para los modelos de asociación transformaremos el dataset original, de tal manera que podamos observar la relación entre las distintas variables de diagnósticos médicos del dataset.
Definiremos unos objetivos en función del modelo que vamos a utilizar.
A través de los datos del dataset tendremos que llevar a cabo distintos métodos de clasificación donde, y con la mayor eficiencia posible, tendremos que determinar que pacientes padecerán una enfermedad cardíaca o no.
Por otro lado, realizaremos distintos métodos de clustering para estudiar los posibles grupos existentes de aquellos pacientes y gracias a las variables que no han sido previamente categorizadas.
Respecto a la asociación, estudiaremos para los pacientes objetivo, cuales son las constantes médicas que suelen estar asociadas.
Este proyecto lo dividiremos en dos partes diferenciadas, ya que para poder llevar a cabo los objetivos anteriores, en primer lugar deberemos conocer nuestros datos, realizar las transformaciones necesarias, etc. Obtendremos un dataframe distinto para cada modelo a aplicar.
En la segunda parte (PRAC_2), desarrollaremos los modelos anteriores y obtendremos las conclusiones finales del modelo.
Instalamos y cargamos las librerías ggplot2 y dplry
# https://cran.r-project.org/web/packages/ggplot2/index.html
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')
# https://cran.r-project.org/web/packages/dplyr/index.html
if (!require('dplyr')) install.packages('dplyr'); library('dplyr')Instalamos y cargamos la librería knitr
if (!require('knitr')) install.packages('knitr'); library('knitr')Instalamos y cargamos la librería factorextra
if (!require('factoextra')) install.packages('factoextra'); library('factoextra')Instalamos y cargamos la librería tidyr y readr
if (!require('tidyr')) install.packages('tidyr'); library('tidyr')## Loading required package: tidyrif (!require('readr')) install.packages('readr'); library('readr')## Loading required package: readrEste dataset dispone de 4 dataframes sobre el diagnostico de enfermedades cardíacas, los datos se recopilaron de distintas instituciones y por ello proceden de cuatro fuentes diferentes:
1. Cleveland Clinic Foundation (cleveland.data)
2. Instituto Húngaro de Cardiología, Budapest (hungarian.data)
3. V.A. Medical Center, Long Beach, CA (long-beach-va.data)
4. Hospital Universitario, Zúrich, Suiza (switzerland.data)
Cada dataframe contiene el mismo formato, utilizando las mismas variables o atributos, 14, aunque cabe resaltar que el dataset original disponía de 76 variables, nosotros en este estudio utilizaremos los 14.
Este dataset y por ello este proyecto ha sido posible gracias a los siguientes investigadores:
1. Instituto Húngaro de Cardiología, Budapest: Andras Janosi, M.D.
2. Hospital Universitario, Zúrich, Suiza: William Steinbrunn, M.D.
3. Hospital Universitario, Basilea, Suiza: Matthias Pfisterer, M.D.
4. V.A. Medical Center, Long Beach y Cleveland Clinic Foundation: Robert Detrano, M.D., Ph.D.El dataset hace referencia a las siguientes observaciones para cada institución: Cleveland: 303 Hungarian: 294 Switzerland: 123 Long Beach VA: 200
En total, 920 observaciones.
Entre ellas encontramos:
age: edad, variable continua.
Gender: Sexo de cada individuo. 1 = Hombre, 0 = Mujer, variable categórica binaria.
cp: Tipo de dolo de pecho (chest pain), variable categórica, valores numéricos del 1 al 4, donde:
1: Agina típica
2: Angina atípica
3: Dolor no anginoso
4: Asintomaticotrestbps: presión arterial en reposo (en mm Hg al ingreso al hospital), variable numérica continua.
chol: Colesterol sérico en mg / dl, variable continua. Los niveles de colesterol total se pueden clasificar según criterios médicos de la siguiente manera: Deseables: < 200 mg/dl. Límite alto: 200-239 mg/dl. Alto: ≥ 240 mg/dl.
FBS: Azúcar en sangre en ayunas superior a 120 mg / dl, variable categórica binaria, 1 = si, 0 = no.
RestECG: Resultados electrocardiográficos en reposo. Variable categórica numérica con valores 0, 1 y 2.
0: normal
1: tener anomalía de la onda ST-T (inversiones de la onda T y / o ST elevación o depresión de> 0,05 mV).
2: muestra hipertrofia ventricular izquierda probable o definitivaThalach: frecuencia cardíaca máxima alcanzada, valor numérico continuo.
Exang: Angina inducida por ejercicio si o no, variable categórica binaria, 1 = si, 0 = no.
Oldpeak: Depresión del ST inducida por el ejercicio en relación con el reposo. Variable continua numérica float.
Slope: Pendiente del pico de ejercicio segmento ST, variable categórica numérica con valores 1, 2 y 3.
1: pendiente ascendente
2: plano
3: pendiente descendenteCA: número de vasos principales coloreados por fluoroscopia, variable categórica numérica de 0 a 3.
Thal: indica la talasemia, variable categórica numérica con tres valores: 3,6 y 7.
3: normal
6: Defecto fijo
7: Defecto reversible.Diagnosis: Variable objetivo. Muestra si la persona padece una enfermedad cardíaca o no:
0: ausencia
1, 2, 3, 4: presente.Al tratarse de variables tan especificas, en el ámbito de la cardiología realizaremos un estudio más exhaustivo de cada una para poder comprenderlas mejor.
Age : La edad es el factor de riesgo más importante en el desarrollo de enfermedades cardiovasculares o cardíacas, según la SEC, con aproximadamente un triple de riesgo con cada década de vida. Las vetas de grasa coronarias pueden comenzar a formarse en la adolescencia. Se estima que el 82 por ciento de las personas que mueren de enfermedad coronaria tienen 65 años o más. Simultáneamente, el riesgo de accidente cerebrovascular se duplica cada década después de los 55 años.
Sex: Los hombres tienen un mayor riesgo de enfermedad cardíaca que las mujeres premenopáusicas. Una vez pasada la menopausia, se ha argumentado que el riesgo de una mujer es similar al del hombre, aunque datos más recientes de la OMS y la ONU lo niegan. Si una mujer tiene diabetes, tiene más probabilidades de desarrollar una enfermedad cardíaca que un hombre con diabetes.
Cp (Angina o dolor de pecho): La angina es un dolor o malestar en el pecho que se produce cuando el músculo cardíaco no recibe suficiente sangre rica en oxígeno.
Trestbps: la presión arterial alta puede dañar las arterias que alimentan el corazón. La presión arterial alta puede ocurrir debido a distintas afecciones, como la obesidad, el colesterol alto o la diabetes.
Chol: Un nivel alto de colesterol unido a lipoproteínas de baja densidad (LDL) (el colesterol “malo”) ayuda a que las arterias se estrechen. Un nivel alto de triglicéridos, un tipo de grasa en la sangre, también aumenta el riesgo de un ataque cardíaco. Sin embargo, un nivel alto de colesterol de lipoproteínas de alta densidad (HDL) (el colesterol “bueno”) reduce el riesgo de sufrir un ataque cardíaco.
FBS: No producir suficiente hormona secretada por el páncreas (insulina) o no responder a la insulina de manera adecuada hace que los niveles de azúcar en sangre del cuerpo aumenten, lo que aumenta el riesgo de un ataque cardíaco.
RestECG: Para las personas con bajo riesgo de enfermedad cardiovascular, la USPSTF concluye con certeza moderada que los daños potenciales de la detección con ECG en reposo o de ejercicio igualan o superan los beneficios potenciales. Para las personas con riesgo intermedio a alto, la evidencia actual es insuficiente para evaluar el balance de beneficios y daños del cribado.
Thalach: El aumento del riesgo cardiovascular, asociado con la aceleración de la frecuencia cardíaca, fue comparable al aumento del riesgo observado con la presión arterial alta. Se ha demostrado que un aumento de la frecuencia cardíaca de 10 latidos por minuto se asoció con un aumento del riesgo de muerte cardíaca en al menos un 20%, y este aumento del riesgo es similar al observado con un aumento de la sangre sistólica. presión por 10 mm Hg.
Oldpeak : El infarto de miocardio con elevación del segmento ST es una emergencia médica producida por la formación de un trombo sobre una placa rota de aterosclerosis que ocluye la circulación coronaria del músculo cardíaco.
Slope: Pendiente del pico de ejercicio segmento ST,
CA: Un vaso sanguíneo es una estructura hueca y tubular que conduce la sangre impulsada por la acción del corazón, cuya función principal es transportar nutrientes, oxígeno y desechos del cuerpo. Existen cinco tipos de vasos sanguíneos; los cuales se clasifican en arterias, arteriolas, capilares, vénulas y venas (ordenados por el recorrido que realiza la sangre desde que sale del corazón hasta que retorna al mismo). Los vasos sanguíneos forman parte del aparato circulatorio. Los tres tipos principales de vasos sanguíneos son: Arterias, venas y capilares. Los vasos sanguíneos transportan sangre por todo el cuerpo. Las arterias transportan sangre desde el corazón. Las venas llevan la sangre de regreso al corazón. Los capilares rodean a las células y a los tejidos del cuerpo para aportar y absorber oxígeno, nutrientes y otras sustancias. Los capilares también conectan las ramas de las arterias y las ramas de las venas.
Thal: La talasemia es un trastorno de la sangre hereditario (es decir, se pasa de los padres a los hijos a través de los genes) que ocurre cuando el cuerpo no produce la cantidad suficiente de una proteína llamada hemoglobina, una parte importante de los glóbulos rojos.
Creamos un directorio data, en caso de que no exista, para guardar los ficheros que vamos a utilizar como fuente de datos y otro para guardar los gráficos.
if(!file.exists("data")) {
dir.create("data")
}
if(!file.exists("gráficos")) {
dir.create("gráficos")
}Guardamos los ficheros en el directorio data:
URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data"
download.file(URL,destfile = "data/cleveland.data")
URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.hungarian.data"
download.file(URL,destfile = "data/hungarian.data")
URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.switzerland.data"
download.file(URL,destfile = "data/switzerland.data")
URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.va.data"
download.file(URL,destfile = "data/va.data")Cargamos todos los datos, los agrupamos y añadimos los nombres de los atributos
cleve <- read.table('data/cleveland.data', sep = ',', na.strings = "?", fill = F, strip.white = T)
hung <- read.table('data/hungarian.data', sep = ',',na.strings = "?", fill = F, strip.white = T)
swiss <- read.table('data/switzerland.data', sep = ',', na.strings = "?", fill = F, strip.white = T,)
va <- read.table('data/va.data', sep = ',', fill = F, na.strings = "?", strip.white = T)
df <- rbind(cleve, hung, swiss, va)
names(df) <- c('Age', 'Gender', 'CP', 'Trestbps', 'Chol', 'FBS', 'RestECG',
'Thalach', 'Exang', 'Oldpeak', 'Slope', 'CA', 'Thal', 'Diagnosis')Verificamos la estructura del juego de datos principal:
str(df)## 'data.frame': 920 obs. of 14 variables:
## $ Age : num 63 67 67 37 41 56 62 57 63 53 ...
## $ Gender : num 1 1 1 1 0 1 0 0 1 1 ...
## $ CP : num 1 4 4 3 2 2 4 4 4 4 ...
## $ Trestbps : num 145 160 120 130 130 120 140 120 130 140 ...
## $ Chol : num 233 286 229 250 204 236 268 354 254 203 ...
## $ FBS : num 1 0 0 0 0 0 0 0 0 1 ...
## $ RestECG : num 2 2 2 0 2 0 2 0 2 2 ...
## $ Thalach : num 150 108 129 187 172 178 160 163 147 155 ...
## $ Exang : num 0 1 1 0 0 0 0 1 0 1 ...
## $ Oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ Slope : num 3 2 2 3 1 1 3 1 2 3 ...
## $ CA : num 0 3 2 0 0 0 2 0 1 0 ...
## $ Thal : num 6 3 7 3 3 3 3 3 7 7 ...
## $ Diagnosis: int 0 2 1 0 0 0 3 0 2 1 ...Estadísticas de valores vacíos:
colSums(is.na(df))## Age Gender CP Trestbps Chol FBS RestECG Thalach
## 0 0 0 59 30 90 2 55
## Exang Oldpeak Slope CA Thal Diagnosis
## 55 62 309 611 486 0Como observamos, encontramos un número elevado de valores nulos en el dataset, por lo cual deberemos estudiar como vamos a tratar estos registros:
Para aquellas variables donde el número de valores nulos es menor que 100 vamos a eliminar dichos registros ya que representarán menos del 10% de los registros del dataframe:
df <- df[!is.na(df$Trestbps),]
df <- df[!is.na(df$Chol),]
df <- df[!is.na(df$FBS),]
df <- df[!is.na(df$RestECG),]
df <- df[!is.na(df$Thalach),]
df <- df[!is.na(df$Exang),]
df <- df[!is.na(df$Oldpeak),]Observamos el resumen del dataframe resultante:
str(df)## 'data.frame': 740 obs. of 14 variables:
## $ Age : num 63 67 67 37 41 56 62 57 63 53 ...
## $ Gender : num 1 1 1 1 0 1 0 0 1 1 ...
## $ CP : num 1 4 4 3 2 2 4 4 4 4 ...
## $ Trestbps : num 145 160 120 130 130 120 140 120 130 140 ...
## $ Chol : num 233 286 229 250 204 236 268 354 254 203 ...
## $ FBS : num 1 0 0 0 0 0 0 0 0 1 ...
## $ RestECG : num 2 2 2 0 2 0 2 0 2 2 ...
## $ Thalach : num 150 108 129 187 172 178 160 163 147 155 ...
## $ Exang : num 0 1 1 0 0 0 0 1 0 1 ...
## $ Oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ Slope : num 3 2 2 3 1 1 3 1 2 3 ...
## $ CA : num 0 3 2 0 0 0 2 0 1 0 ...
## $ Thal : num 6 3 7 3 3 3 3 3 7 7 ...
## $ Diagnosis: int 0 2 1 0 0 0 3 0 2 1 ...Tenemos 740 observaciones con 14 variables.
Para variable slope observamos 229 valores nulos, siendo los valores usados en esta variable, 1: ,2: y 4:
Para variable CA observamos 463 valores nulos, siendo los valores usados en esta variable, 0 a 3, número de vasos principales coloreados por fluoroscopia.
Y para variable Thal observamos 366 valores nulos, siendo los valores usados en esta variable, 3: normal, 6: Defecto fijo, 7: Defecto reversible, indicadores de la talasemia.
Como podemos ver el número de registros con valores nulos es muy amplio, no podremos eliminarlos ya que en algunos casos alcanza hasta casi un 50% de las muestras del modelo. Por lo tanto, para todas estas variables vamos a asociar el valor numérico “10” que indicará que es un dato desconocido. 10 = Desconocido
df[is.na(df)] <- 10Observamos el resumen del dataframe resultante:
summary(df) ## Age Gender CP Trestbps
## Min. :28.0 Min. :0.0000 Min. :1.000 Min. : 0.0
## 1st Qu.:46.0 1st Qu.:1.0000 1st Qu.:2.000 1st Qu.:120.0
## Median :54.0 Median :1.0000 Median :4.000 Median :130.0
## Mean :53.1 Mean :0.7649 Mean :3.227 Mean :132.8
## 3rd Qu.:60.0 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:140.0
## Max. :77.0 Max. :1.0000 Max. :4.000 Max. :200.0
## Chol FBS RestECG Thalach Exang
## Min. : 0.0 Min. :0.00 Min. :0.0000 Min. : 60.0 Min. :0.0
## 1st Qu.:197.0 1st Qu.:0.00 1st Qu.:0.0000 1st Qu.:120.0 1st Qu.:0.0
## Median :231.0 Median :0.00 Median :0.0000 Median :140.0 Median :0.0
## Mean :220.1 Mean :0.15 Mean :0.6351 Mean :138.7 Mean :0.4
## 3rd Qu.:271.0 3rd Qu.:0.00 3rd Qu.:1.0000 3rd Qu.:159.2 3rd Qu.:1.0
## Max. :603.0 Max. :1.00 Max. :2.0000 Max. :202.0 Max. :1.0
## Oldpeak Slope CA Thal
## Min. :-1.0000 Min. : 1.000 Min. : 0.000 Min. : 3.000
## 1st Qu.: 0.0000 1st Qu.: 2.000 1st Qu.: 1.000 1st Qu.: 3.000
## Median : 0.5000 Median : 2.000 Median :10.000 Median : 7.000
## Mean : 0.8943 Mean : 4.091 Mean : 6.177 Mean : 7.315
## 3rd Qu.: 1.5000 3rd Qu.:10.000 3rd Qu.:10.000 3rd Qu.:10.000
## Max. : 6.2000 Max. :10.000 Max. :10.000 Max. :10.000
## Diagnosis
## Min. :0.0000
## 1st Qu.:0.0000
## Median :1.0000
## Mean :0.9243
## 3rd Qu.:1.0000
## Max. :4.0000En este momento ya hemos resuelto la problemática de los valores nulos, pero deberemos recordar durante todo el modelo que para las variables: CA, Thal y Slope el valor 10 corresponde a los valores desconocidos.
Vamos a utilizar el gráfico boxplot para estudiar la existencia de outliers en el modelo y la distribución de dichas variables, en concreto en las variables continuas.
plAge <- ggplot(df, aes(x=factor(Diagnosis),y=Age))
plAge + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw()
means <- aggregate(Trestbps ~ Diagnosis,df,mean)
plTrestbps <- ggplot(df, aes(x=factor(Diagnosis),y=Trestbps))
plTrestbps + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw() +
geom_text(data = means, aes(label = Trestbps, y = Trestbps + 0.08), size = 2.5)
plChol <- ggplot(df, aes(x=factor(Diagnosis),y=Chol))
plChol + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw()
plThalach <- ggplot(df, aes(x=factor(Diagnosis),y=Thalach))
plThalach + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw()
plThalach <- ggplot(df, aes(x=factor(Diagnosis),y=Oldpeak))
plThalach + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw()
Debemos eliminar aquel registro, outlier que tiene la presión arterial igual a a cero
df <- df[!(df$Trestbps == 0), ] means <- aggregate(Trestbps ~ Diagnosis,df,mean)
plTrestbps <- ggplot(df, aes(x=factor(Diagnosis),y=Trestbps))
plTrestbps + geom_boxplot(aes(fill=factor(Diagnosis))) + theme_bw() +
geom_text(data = means, aes(label = Trestbps, y = Trestbps + 0.08), size = 2.5)
Observamos que la media ha aumentado.
Este es el momento adecuado para crear el dataframe que utilizaremos en el clustering, en la parte dos del proyecto decidiremos que variables utilizar.
dfk <- df
dfk$Age <- NULL
dfk$Gender <- NULL
dfk$FCMax <- NULL
dfk$CP <- NULL
dfk$FBS <- NULL
dfk$Slope <- NULL
dfk$Exang <- NULL
dfk$RestECG <- NULL
dfk$CA <- NULL
dfk$Thal <- NULL
dfk$Diagnosis <- NULL head(dfk)## Trestbps Chol Thalach Oldpeak
## 1 145 233 150 2.3
## 2 160 286 108 1.5
## 3 120 229 129 2.6
## 4 130 250 187 3.5
## 5 130 204 172 1.4
## 6 120 236 178 0.8Antes de continuar vamos a modificar la variable objetivo y simplificarla, la convertiremos a una variable categórica binaria donde, 0: No padece ninguna enfermedad cardíaca y 1: Si que padece, que agrupara los valores del 1 al 4.
df$Diagnosis <- gsub(0, 0, df$Diagnosis)
df$Diagnosis <- gsub(1, 1, df$Diagnosis)
df$Diagnosis <- gsub(2, 1, df$Diagnosis)
df$Diagnosis <- gsub(3, 1, df$Diagnosis)
df$Diagnosis <- gsub(4, 1, df$Diagnosis)
df$Diagnosis<- as.factor(df$Diagnosis)A esta variable la llamaremos FCMax y será la frecuencia máxima , que calcularemos con la siguiente formula (FCMax = 220 - Edad). Corresponde a la fórmula de Fox y Haskell, una manera estándar y con cierto margen de error, pero puede sernos útil para comparar la nueva variable con la variable Thalach .
Comparando las dos variables, crearemos una nueva variable que indica si la Frecuencia máxima de cada persona ha sido superada o no. Esta nueva variable se denominará FCmax_Thalach y será binaria: 1 : Ha sido superada, 0: No ha sido superada.
Creamos las nuevas variables:
df$FCMax <- 220 - df$Age
## Transformación de variables:df$FCMax_Thalac <- df$FCMax - df$Thalach
## Transformación de variables:Sustituimos los valores por 0 y 1
df$FCMax_Thalac <- with(df, ifelse(FCMax_Thalac <= 0, 0,1)) Borramos la variables que no necesitamos:
df$FCMax <- NULL
df$Thalach<- NULLhead(df)## Age Gender CP Trestbps Chol FBS RestECG Exang Oldpeak Slope CA Thal Diagnosis
## 1 63 1 1 145 233 1 2 0 2.3 3 0 6 0
## 2 67 1 4 160 286 0 2 1 1.5 2 3 3 1
## 3 67 1 4 120 229 0 2 1 2.6 2 2 7 1
## 4 37 1 3 130 250 0 0 0 3.5 3 0 3 0
## 5 41 0 2 130 204 0 2 0 1.4 1 0 3 0
## 6 56 1 2 120 236 0 0 0 0.8 1 0 3 0
## FCMax_Thalac
## 1 1
## 2 1
## 3 1
## 4 0
## 5 1
## 6 0Vamos a crear nuevas variables transformando las variables continuas en discretas, lo haremos estudiando su distribución, comparándolas con la variable objetivo y estudiando los criterios médicos.
Los criterios médicos son los siguientes:
Variable Age, > 65 alto riesgo, >55 riesgo medio, <55 riesgo bajo.
La variable chol : Se considera deseable un nivel de colesterol menor de 200 mg/dl, alto entre 200 y 239 y muy alto mayor de 240 mg/dl.
Estudiamos la distribución de las distintas variables, haciendo énfasis en aquellas variables numéricas continuas, que procederemos a discretizar para llevar a cabo un estudio más exhaustivo.
ggplot(data=df,aes(x=Age))+geom_bar() + ggtitle("Número de registros por edad")
ggplot(data=df,aes(x=Chol))+geom_bar() + ggtitle("Número de registros por coresterol")
ggplot(data=df,aes(x=Trestbps))+geom_bar() + ggtitle("Número de registros por presión arterial")
ggplot(data=df,aes(x=Oldpeak))+geom_bar() + ggtitle("Número de registros por presión arterial")
Comenzamos por la variable edad:
summary(df[,"Age"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.00 46.00 54.00 53.09 60.00 77.00Teniendo en cuenta los criterios médicos y la distribución de los registros hemos decidido discretizarlos en los distintos grupos:
0-39(Riesgo muy bajo)
40-54(Riesgo bajo)
55-64(Riesgo medio)
65-100(Riesgo alto)df["edad_discretizada"] <- cut(df$Age, breaks = c(0,40,55,65,100), labels = c("0-39(Riesgo muy bajo)", "40-54(Riesgo bajo)", "55-64(Riesgo medio)", "65-100(Riesgo alto)"))Variable coresterol, en este caso vamos a discretizar la variable únicamente siguiendo criterios médicos, ya que han sido rigurosamente establecidos.
summary(df[,"Chol"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 197.0 231.0 220.4 271.0 603.0Se considera deseable un nivel de colesterol menor de 200 mg/dl, alto entre 200 y 239 y muy alto mayor de 240 mg/dl
df["chol_discretizado"] <- cut(df$Chol, breaks = c(-1,200,240,1000), labels = c("0-199(Riesgo bajo)", "200-239(Riesgo medio)", "240-x(Riesgo alto)"))Variable Trestbps, Presión arterial:
Para esta variable también tendremos en cuenta las categorías establecidas por los cardiólogos:
Normal, menos de 120
Elevada: 120-129
Presión arterial alta nivel 1: 130-139
Presión arterial alta nivel 2: 140-179
Crisis de hipertensión: mayor de 180df["trestbps_discretizado"] <- cut(df$Trestbps, breaks = c(0,120,130,140,180,500), labels = c("0-119(Normal)", "120-129(Elevada)", "130-139(Alta nivel 1)", "140-179(Alta nivel 2)", "180-x(Crisis de hipertesión)"))variable Oldpeak:
df["oldpeak_discretizado"] <- cut(df$Oldpeak, breaks = c(-2,0,2,4,6,8,10), labels = c(0, 1, 2, 3,4,5))
df$oldpeak_discretizado <- as.numeric(as.character(df$oldpeak_discretizado))Observamos el dataframe resultante:
head(df)## Age Gender CP Trestbps Chol FBS RestECG Exang Oldpeak Slope CA Thal Diagnosis
## 1 63 1 1 145 233 1 2 0 2.3 3 0 6 0
## 2 67 1 4 160 286 0 2 1 1.5 2 3 3 1
## 3 67 1 4 120 229 0 2 1 2.6 2 2 7 1
## 4 37 1 3 130 250 0 0 0 3.5 3 0 3 0
## 5 41 0 2 130 204 0 2 0 1.4 1 0 3 0
## 6 56 1 2 120 236 0 0 0 0.8 1 0 3 0
## FCMax_Thalac edad_discretizada chol_discretizado
## 1 1 55-64(Riesgo medio) 200-239(Riesgo medio)
## 2 1 65-100(Riesgo alto) 240-x(Riesgo alto)
## 3 1 65-100(Riesgo alto) 200-239(Riesgo medio)
## 4 0 0-39(Riesgo muy bajo) 240-x(Riesgo alto)
## 5 1 40-54(Riesgo bajo) 200-239(Riesgo medio)
## 6 0 55-64(Riesgo medio) 200-239(Riesgo medio)
## trestbps_discretizado oldpeak_discretizado
## 1 140-179(Alta nivel 2) 2
## 2 140-179(Alta nivel 2) 1
## 3 0-119(Normal) 2
## 4 120-129(Elevada) 2
## 5 120-129(Elevada) 1
## 6 0-119(Normal) 1Una vez que tenemos las variables discretizadas, vamos transformar sus registros a valores numéricos, de esta manera nos serán útiles para realizar los modelos de clasificación o el modelo PCA, que realizaremos a continuación.
Representamos gráficamente las nuevas variables discretizadas y el resto de variables:
Variables discretizadas:
ggplot(data=df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar() + ggtitle("Tipo de trastorno por edad(G1)") + labs(x="Edad", colour= "Diagnóstico") 
ggplot(data=df,aes(x=chol_discretizado,fill=Diagnosis))+geom_bar() + ggtitle("Tipo de trastorno por coresterol(G2)") + labs(x="Nivel de Coresterol", colour= "Diagnóstico")
ggplot(data=df,aes(x=trestbps_discretizado,fill=Diagnosis))+geom_bar() + ggtitle("Tipo de trastorno por presión arterial(G3)") + labs(x="Presión Arterial", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
ggplot(data=df,aes(x=oldpeak_discretizado,fill=Diagnosis))+geom_bar() + ggtitle("Tipo de trastorno por oldpeak (G4)") + labs(x="Oldpeak", colour= "Diagnóstico")
G1 - La variable edad influye en el riesgo de padecer enfermedades cardíacas o no, cuanto mayor es la edad del paciente mayor es la probabilidad de sufrir una enfermedad cardíaca.
G2 - Cuanto mayor es el nivel de colesterol más riesgo existe de sufrir estas enfermedades aunque no es una variable decisiva.
G3 - La presión arterial es un factor de riego importante para prevenir estas enfermedades, siendo mayor de 140 un factor crítico.
G4 - Sin duda, gráficamente podemos observar que la variable discretizada más relevante es oldpeak, cuando la depresión del ST es mayor que cero el riesgo es crítico
En este punto, vamos a realizar el análisis PCA y volveremos al análisis gráfico una vez estudiadas las variables más importantes del modelo.
Vamos a utilizar la variable edad discretizada ya que es una variable con peso en el análisis pca, además de una variable que facilitará nuestro análisis.
Vamos a realizar el análisis gráfico de las variables más relevantes en referencia al análisis PCA y de la variable sexo en relación con la edad,
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~Gender)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y Sexo (G1)") + labs(x="Gender", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
G1 - Los hombres tienen más probabilidad de sufrir enfermedades cardíacas y esta aumenta conforme su edad aumenta.
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~CP)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y Tipo de dolor de pecho (G2)") + labs(x="Tipo de dolor", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
G2 - Aquellas personas asintomáticas (4), seguidas de un dolor no anginoso son las más vulnerables. No tener síntomas de este tipo, no garantiza no sufrir una enfermedad cardíaca.
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~Slope)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y pendiente del pico del segmento (G3)") + labs(x="Edad", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~oldpeak_discretizado)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y depresión del oldpeak (G4)") + labs(x="Edad", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~Exang)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y angina inducida por el ejercicio 1:si 0:no (G5)") + labs(x="Edad", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
ggplot(data = df,aes(x=edad_discretizada,fill=Diagnosis))+geom_bar(position="fill")+facet_wrap(~trestbps_discretizado)+ theme(text = element_text(size=8)) + ggtitle("Diagnóstico por Edad y Presión arterial en reposo (G6)") + labs(x="Edad", y="", colour= "Diagnóstico") + theme(text = element_text(size=12), axis.text.x = element_text(angle=90, hjust=1))
G3 - Observando la relación entre Slope y la edad, concluimos que una pendiente del pico del segmento ST, plana es sinónimo de enfermedad cardíaca para todos los grupos de edad y una pendiente descendente también confirma un riesgo de enfermedad cradíaca, especialmente para las edades más avanzadas.
G4 - La variable oldpeak, es determinante, para valores mayores que uno el riesgo para todas las edades es mayor al 50%, si la depresión es mayor o igual que 2 estamos hablando de un riesgo mayor que el 70%.
G5 - Si la angina ha sido enducida por el ejercicio, el riesgo de padecer una enfermedad cardíaca es muy alta para todas las edades, incrementandose cuando esta aumenta.
G6 - La variable presión arterial en reposo es importante a la hora de padecer enfermedades cardíacas en edades superiores a los 50 años y conforme la edad aumenta y la presión arterial también lo hace este riesgo es mayor.
En este momento vamos a guardar dos dataframes diferenciados para distintos tipos de análisis.
df2 <- dfVamos a sustituir las variables categóricas con formato string por variables numéricas y eliminar todas las variables numéricas:
Edad:
0: 0-39(Riesgo muy bajo)
1: 40-54(Riesgo bajo)
2: 55-64(Riesgo medio)
3: 65-100(Riesgo alto)Colesterol:
0: 0-199(Riesgo bajo)
1: 200-239(Riesgo medio)
2: 240-x(Riesgo alto)Presión arterial
0: 0-119(Normal)
1: 120-129(Elevada)
2: 130-139(Alta nivel 1)
3: 140-179(Alta nivel 2)
4: 180-x(Crisis de hipertesión)df2$edad_discretizada <- with(df2, ifelse(edad_discretizada == "0-39(Riesgo muy bajo)", 0, ifelse(edad_discretizada == "40-54(Riesgo bajo)", 1, ifelse(edad_discretizada == "55-64(Riesgo medio)", 2, 3))))df2$chol_discretizado <- with(df2, ifelse(chol_discretizado == "0-199(Riesgo bajo)", 0, ifelse(chol_discretizado == "200-239(Riesgo medio)", 1, 2)))df2$trestbps_discretizado <- with(df2, ifelse(trestbps_discretizado == "0-119(Normal)", 0, ifelse(trestbps_discretizado == "120-129(Elevada)", 1, ifelse(trestbps_discretizado == "130-139(Alta nivel 1)",2,ifelse(trestbps_discretizado == "140-179(Alta nivel 2)", 3, 4)))))Para reducir la dimensionalidad del dataframe vamos a utilizar el modelo PCA, Principal Component Analysis. Además utilizaremos este modelo para estudiar aquellas variables que tienen más contribución en el modelo
En nuestro caso aplicaremos el método SVD, Single Value Decoposition
PCA es una técnica de extracción de características donde combinamos las entradas de una manera específica y podemos eliminar algunas de las variables “menos importantes” manteniendo la parte más importante de todas las variables. Como valor añadido, aplicando PCA conseguiremos que todas las nuevas variables sean independientes una de otra.
Para realizar el modelo PCA deberemos llevar a cabo los siguientes pasos:
Estandarizar los datos de entrada
Obtener los autovectores y autovalores de la matriz de covarianza
Ordenar los autovalores de mayor a menor y elegir los “k” autovectores que se correspondan con los autovectores “k” más grandes (donde “k” es el número de dimensiones del nuevo subespacio de características).
Construir la matriz de proyección W con los “k” autovectores seleccionados. Transformamos el dataset original “X estandarizado” vía W para obtener las nuevas características k-dimensionales.
Todas nuestras variables categóricas son numéricas por lo que podremos usar el algoritmo de PCA ya que solo puede aplicarse a datos numéricos. Si los datos contienen variables categóricas string, deberian ser convertidas a numéricas, ya hemos realizado este proceso anteriormente.
Utilizaremos la función prcomps ya que junto la función PCA() usan el modelo SVD y realizan todos los pasos descritos anteriormente de forma directa.
head(df2)## Age Gender CP Trestbps Chol FBS RestECG Exang Oldpeak Slope CA Thal Diagnosis
## 1 63 1 1 145 233 1 2 0 2.3 3 0 6 0
## 2 67 1 4 160 286 0 2 1 1.5 2 3 3 1
## 3 67 1 4 120 229 0 2 1 2.6 2 2 7 1
## 4 37 1 3 130 250 0 0 0 3.5 3 0 3 0
## 5 41 0 2 130 204 0 2 0 1.4 1 0 3 0
## 6 56 1 2 120 236 0 0 0 0.8 1 0 3 0
## FCMax_Thalac edad_discretizada chol_discretizado trestbps_discretizado
## 1 1 2 1 3
## 2 1 3 2 3
## 3 1 3 1 0
## 4 0 0 2 1
## 5 1 1 1 1
## 6 0 2 1 0
## oldpeak_discretizado
## 1 2
## 2 1
## 3 2
## 4 2
## 5 1
## 6 1Convertimos la variable diagnosis a un formato numérico
df2$Diagnosis <- as.numeric(as.character(df2$Diagnosis))Función del PCA SVD
pca <- prcomp(df2, scale = TRUE)Gráficos PCA:
fviz_screeplot(pca, addlabels = TRUE)
fviz_contrib(pca, addlabels = TRUE, choice = "var", axes = 1, top = 20)
La línea roja indica el valor medio de contribución. Podemos decir que aquellas variables superiores a esta linea son importantes a la hora de contribuir a este componente.
fviz_pca_ind(pca,
col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)## Warning: ggrepel: 671 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Representa las observaciones mediante sus proyecciones
fviz_pca_var(pca,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
Representa las variables sobre las dos primeras componentes principales. La correlación entre una variable y una componente principal se utiliza como la coordenada de dicha variable sobre la componente principal. De esta manera podemos obtener un gráfico de correlación de variables
fviz_pca_biplot(pca, repel = TRUE,
col.var = "black",
col.ind = "red"
)## Warning: ggrepel: 667 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Representación conjunta de variables y observaciones
Eliminaremos las siguientes variables del dataframe para reducir su dimensionalidad:
FBS , CA,Thal, Chol
Guardamos el dataframe para pode utilizarlo para crear el dataframe que usaremos en el modelo de asociación
df4 <-df2df2$FBS<- NULL
df2$CA<- NULL
df2$Chol<- NULL
df2$Thal<- NULLObservamos el dataframe resultante:
head(df2)## Age Gender CP Trestbps RestECG Exang Oldpeak Slope Diagnosis FCMax_Thalac
## 1 63 1 1 145 2 0 2.3 3 0 1
## 2 67 1 4 160 2 1 1.5 2 1 1
## 3 67 1 4 120 2 1 2.6 2 1 1
## 4 37 1 3 130 0 0 3.5 3 0 0
## 5 41 0 2 130 2 0 1.4 1 0 1
## 6 56 1 2 120 0 0 0.8 1 0 0
## edad_discretizada chol_discretizado trestbps_discretizado
## 1 2 1 3
## 2 3 2 3
## 3 3 1 0
## 4 0 2 1
## 5 1 1 1
## 6 2 1 0
## oldpeak_discretizado
## 1 2
## 2 1
## 3 2
## 4 2
## 5 1
## 6 1Vamos a crear dos dataframes diferentes para los modelos de clasificación, los cuales podremos comparar entre sí:
DataFrame_clasificación_1:
dfc <-df2df3 <- df2Eliminamos las variables que han sido discretizadas
df3$Age<- NULL
df3$Trestbps<- NULL
df3$Oldpeak<- NULLDataFrame_clasificación_2:
dfc1 <- df3head(dfc1)## Gender CP RestECG Exang Slope Diagnosis FCMax_Thalac edad_discretizada
## 1 1 1 2 0 3 0 1 2
## 2 1 4 2 1 2 1 1 3
## 3 1 4 2 1 2 1 1 3
## 4 1 3 0 0 3 0 0 0
## 5 0 2 2 0 1 0 1 1
## 6 1 2 0 0 1 0 0 2
## chol_discretizado trestbps_discretizado oldpeak_discretizado
## 1 1 3 2
## 2 2 3 1
## 3 1 0 2
## 4 2 1 2
## 5 1 1 1
## 6 1 0 1Para el dataframe 3 hemos eliminado todas las variables continuas y de esta manera todas nuestras variables serán categóricas, lo que nos ayudará a realizar lo modelos de clasificación de arboles de decisión de forma eficiente. En la segunda parte compararemos ambos dataframe, dfc y dfc1 observando cuál de ellos lleva a modelos de clasificación más eficientes.
Primero vamos a elegir las variables que vamos a utilizar en el modelo de asociación:
Para este modelo utilizaremos las siguientes variables: CP, RestEG, Slope, Thal, CA y Diagnostico.
En el modelo de asociación, en la segunda parte de este proyecto, estudiaremos la relación de los diferentes diagnósticos médicos procedentes de las variables anteriores
Además añadiremos una nueva variable que llamaremos paciente que identificará a cada individuo por diagnóstico.
Modificaremos las variables y les daremos los siguientes valores:
Partiremos de df4:
head(df4)## Age Gender CP Trestbps Chol FBS RestECG Exang Oldpeak Slope CA Thal Diagnosis
## 1 63 1 1 145 233 1 2 0 2.3 3 0 6 0
## 2 67 1 4 160 286 0 2 1 1.5 2 3 3 1
## 3 67 1 4 120 229 0 2 1 2.6 2 2 7 1
## 4 37 1 3 130 250 0 0 0 3.5 3 0 3 0
## 5 41 0 2 130 204 0 2 0 1.4 1 0 3 0
## 6 56 1 2 120 236 0 0 0 0.8 1 0 3 0
## FCMax_Thalac edad_discretizada chol_discretizado trestbps_discretizado
## 1 1 2 1 3
## 2 1 3 2 3
## 3 1 3 1 0
## 4 0 0 2 1
## 5 1 1 1 1
## 6 0 2 1 0
## oldpeak_discretizado
## 1 2
## 2 1
## 3 2
## 4 2
## 5 1
## 6 1df4$Age<- NULL
df4$Gender<- NULL
df4$Oldpeak<- NULL
df4$Chol<- NULL
df4$FBS<- NULL
df4$Exang<- NULL
df4$Oldpeak<- NULL
df4$FCMax_Thalac<- NULL
df4$edad_discretizada<- NULL
df4$chol_discretizado<- NULL
df4$trestbps_discretizado<- NULL
df4$oldpeak_discretizado<- NULL
df4$Trestbps<- NULLhead(df4)## CP RestECG Slope CA Thal Diagnosis
## 1 1 2 3 0 6 0
## 2 4 2 2 3 3 1
## 3 4 2 2 2 7 1
## 4 3 0 3 0 3 0
## 5 2 2 1 0 3 0
## 6 2 0 1 0 3 0Modificamos las variables numéricas por sus correspondientes diagnósticos médicos
df4$CP <- gsub(1, 'angina_típica', df4$CP)
df4$CP <- gsub(2, 'angina_atípica', df4$CP)
df4$CP <- gsub(3, 'dolor_no_anginoso', df4$CP)
df4$CP <- gsub(4, 'asintomatico', df4$CP)
df4$RestECG <- gsub(0, 'resultados_electro_normal', df4$RestECG)
df4$RestECG <- gsub(1, 'resultados_electro_anomalia', df4$RestECG)
df4$RestECG <- gsub(2, 'resultados_electro_hipertrofia', df4$RestECG)
df4$Slope <- gsub(10, 'desconocido', df4$Slope)
df4$Slope <- gsub(1, 'pendiente_ST_ascendente', df4$Slope)
df4$Slope <- gsub(2, 'pendiente_ST_plano', df4$Slope)
df4$Slope <- gsub(3, 'pendiente_ST_descendente', df4$Slope)
df4$CA <- gsub(10, 'desconocido', df4$CA)
df4$CA <- gsub(0, 'vasos_coloreados_0', df4$CA)
df4$CA <- gsub(1, 'vasos_coloreados_1', df4$CA)
df4$CA <- gsub(2, 'vasos_coloreados_2', df4$CA)
df4$CA <- gsub(3, 'vasos_coloreados_3', df4$CA)
df4$Thal <- gsub(10, 'desconocido', df4$Thal)
df4$Thal <- gsub(3, 'talasemia_normal', df4$Thal)
df4$Thal <- gsub(6, 'talasemia_defecto_fijo', df4$Thal)
df4$Thal <- gsub(7, 'talasemia_defecto_reversible', df4$Thal)
df4$Diagnosis <- gsub(0, 'sin_enfermedad_cardiaca', df4$Diagnosis)
df4$Diagnosis <- gsub(1, 'con_enfermedad_cardiaca', df4$Diagnosis)Creamos una nueva variable que representa el Id del paciente
df4$ID_Paciente <- 1:nrow(df4)head(df4)## CP RestECG Slope
## 1 angina_típica resultados_electro_hipertrofia pendiente_ST_descendente
## 2 asintomatico resultados_electro_hipertrofia pendiente_ST_plano
## 3 asintomatico resultados_electro_hipertrofia pendiente_ST_plano
## 4 dolor_no_anginoso resultados_electro_normal pendiente_ST_descendente
## 5 angina_atípica resultados_electro_hipertrofia pendiente_ST_ascendente
## 6 angina_atípica resultados_electro_normal pendiente_ST_ascendente
## CA Thal Diagnosis
## 1 vasos_coloreados_0 talasemia_defecto_fijo sin_enfermedad_cardiaca
## 2 vasos_coloreados_3 talasemia_normal con_enfermedad_cardiaca
## 3 vasos_coloreados_2 talasemia_defecto_reversible con_enfermedad_cardiaca
## 4 vasos_coloreados_0 talasemia_normal sin_enfermedad_cardiaca
## 5 vasos_coloreados_0 talasemia_normal sin_enfermedad_cardiaca
## 6 vasos_coloreados_0 talasemia_normal sin_enfermedad_cardiaca
## ID_Paciente
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6Para finalizar transformaremos el dataframe, de tal manera que lo podamos usar para el modelo de asociación.
dfa <- pivot_longer(df4, -c(ID_Paciente), values_to = "Value", names_to = "Atributos")dfa$Atributos<- NULL head (dfa)## # A tibble: 6 x 2
## ID_Paciente Value
## <int> <chr>
## 1 1 angina_típica
## 2 1 resultados_electro_hipertrofia
## 3 1 pendiente_ST_descendente
## 4 1 vasos_coloreados_0
## 5 1 talasemia_defecto_fijo
## 6 1 sin_enfermedad_cardiacastr(dfa)## tibble[,2] [4,434 x 2] (S3: tbl_df/tbl/data.frame)
## $ ID_Paciente: int [1:4434] 1 1 1 1 1 1 2 2 2 2 ...
## $ Value : chr [1:4434] "angina_típica" "resultados_electro_hipertrofia" "pendiente_ST_descendente" "vasos_coloreados_0" ...Ya tenemos el dataframe final para el modelo de asociación. Vamos a guardarlo en un csv.
write.csv(dfa,"data/pacientes_diagnosticos.csv", row.names = FALSE)Hemos preparado 3 dataframe distintos, para aplicar cada tipo de modelo dfc y dfc1: Dataframes para realizar los modelos de clasificación. dfk: para realizar los modelos de clustering. dfa: para los modelos de asociación.
Si fuera necesario, en el momento de realizar los modelos, podríamos modificar los dataframes anteriores.
Debido a que este dataset esta diseñado para modelos de clasificación, hemos realizado algunos cambios en los datasets para poder implementar los métodos de clasificación y asociación.
Para el caso de los métodos de clustering, hemos eliminado aquellas variables que por si mismas representan grupos predefinidos, como por ejemplo el sexo y la edad.
Para el caso de los modelos de asociación hemos rehecho completamente el dataframe y hemos creado uno nuevo con dos variables, el ID_paciente y sus diagnósticos médicos , de tal manera que las podremos relacionar en el modelo de asociación.
Abrimos el dataframe generado para los modelos de clustering
head(dfk)## Trestbps Chol Thalach Oldpeak
## 1 145 233 150 2.3
## 2 160 286 108 1.5
## 3 120 229 129 2.6
## 4 130 250 187 3.5
## 5 130 204 172 1.4
## 6 120 236 178 0.8str(dfk)## 'data.frame': 739 obs. of 4 variables:
## $ Trestbps: num 145 160 120 130 130 120 140 120 130 140 ...
## $ Chol : num 233 286 229 250 204 236 268 354 254 203 ...
## $ Thalach : num 150 108 129 187 172 178 160 163 147 155 ...
## $ Oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...Cargamos las librerías necesarias:
if (!require('cluster')) install.packages('cluster'); library('cluster')## Loading required package: clusterCreamos una función para normalizar:
normalizar <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))
}Normalizamos todas las variables: Trestbps, Chol, Thalach y Oldpeak:
dfk$Trestbps_norm<-normalizar(dfk$Trestbps)
dfk$Chol_norm<-normalizar(dfk$Chol)
dfk$Thalach_norm<-normalizar(dfk$Thalach)
dfk$Oldpeak_norm<-normalizar(dfk$Oldpeak)En este caso vamos a llevar a cabo el modelo utilizando el algoritmo K-means.
K-means es un algoritmo de clasificación no supervisada (clusterización) que agrupa objetos en k grupos basándose en sus características. El agrupamiento se realiza minimizando la suma de distancias entre cada objeto y el centroide de su grupo o cluster. Se suele usar la distancia cuadrática.
Puesto que inicialmente, en un problema real, no se conoce cual es el número más idóneo de clústers k, vamos a probar primero con dos (el valor óptimo) y posteriormente de 4 a 8 clústers. Para evaluar la calidad de cada proceso de agregación vamos a usar la silueta media. La silueta de cada muestra evalúa como de bien o mal está clasificada la muestra en el clúster al que ha sido asignada. Para ello se usa una fórmula que tiene en cuenta la distancia a las muestras de su clúster y la distancia a las muestras del clúster vecino más cercano.
A la hora de probar el código que se muestra, es importante tener en cuenta que las muestras se generan de forma aleatoria y también que el algoritmo kmeans tiene una inicialización aleatoria. Por lo tanto, en cada ejecución se obtendrá unos resultados ligeramente diferentes.
Lo primero que hacemos es cargar la librería cluster que contiene las funciones que se necesitan
A continuación vamos a aplicar el algoritmo kmeans de 2 a 8 clústers:
Tras realizar varias pruebas, las variables elegidas para el modelo de clasificación serán las siguientes: Thalach y Trestbps.
Eliminamos también las variables sin normalizar.
dfk$Oldpeak <- NULL
dfk$Thalach<- NULL
dfk$Trestbps<- NULL
dfk$Chol <- NULL
dfk$Oldpeak_norm <- NULL
dfk$Chol_norm <- NULLtotalData <- dfkstr(totalData)## 'data.frame': 739 obs. of 2 variables:
## $ Trestbps_norm: num 0.491 0.63 0.259 0.352 0.352 ...
## $ Thalach_norm : num 0.634 0.338 0.486 0.894 0.789 ...fit2 <- kmeans(totalData, 2)
y_cluster2 <- fit2$cluster
fit3 <- kmeans(totalData, 3)
y_cluster3 <- fit3$cluster
fit4 <- kmeans(totalData, 4)
y_cluster4 <- fit4$cluster
fit5 <- kmeans(totalData, 5)
y_cluster5 <- fit5$cluster
fit6 <- kmeans(totalData, 6)
y_cluster6 <- fit6$cluster
fit7 <- kmeans(totalData, 7)
y_cluster7 <- fit7$cluster
fit8 <- kmeans(totalData, 8)
y_cluster8 <- fit8$clusterLas variables y_cluster2 a y_cluster8 contienen para cada muestra el identificador del clúster a las que han sido asignadas. Por ejemplo, en el caso de los k=2 las muestras se han asignado al clúster 1 o al 2 y para el k=8 de 1 al 8
y_cluster2## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 2 2 2 2 1 2 2 1 2 1 2 2 2 2 2 2 1 1 1 1
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 1 2 2 2 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 2 2 1 2 1 2 2 2 1 2 2 2 1 2 2 2 2 2 2 1
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 2 2 2 2 2
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 1 1
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 2 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 1 1 2 2
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 2 2 2 1 2 2 1 2 2 2 2 2 2 1 1 1 2 2 2 2
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 2 2 2 2 2 2 2 2 2 2 1 1 1 2 1 1 2 1 2 2
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 2 2 2 1 2 2 2 1 2 2 2 1 2 1 1 2 1 2 2 1
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
## 2 1 2 1 2 2 1 1 2 2 2 2 2 2 2 2 2 2 1 2
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## 2 2 2 1 2 2 2 2 1 2 2 1 1 1 2 1 1 2 2 2
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
## 2 2 2 2 1 1 2 1 2 2 2 1 1 2 2 2 2 1 1 2
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
## 2 2 2 2 1 1 2 2 2 2 1 1 1 1 2 2 2 2 2 1
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
## 2 2 1 2 2 2 1 2 2 2 2 2 2 2 2 2 1 1 2 2
## 301 302 303 304 305 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
## 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 322 323 324 325 326 327 328 329 330 332 333 334 336 337 339 340 341 342 343 344
## 2 2 2 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 2 2
## 345 346 347 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
## 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 2 2 2 2
## 366 367 368 370 371 372 373 374 375 377 378 380 381 382 383 384 386 387 388 389
## 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 2 1 1
## 391 392 393 396 397 398 399 400 402 403 404 407 408 409 410 413 414 415 416 417
## 2 2 1 1 1 2 1 1 2 1 2 2 2 2 2 2 1 2 2 1
## 418 419 420 421 422 423 424 425 426 427 429 430 431 432 433 434 436 437 439 440
## 1 1 2 1 2 1 2 2 1 2 2 2 1 1 2 2 2 2 2 2
## 441 442 443 444 445 446 447 449 450 451 452 453 454 455 456 457 459 460 461 462
## 1 2 2 2 2 2 2 1 1 2 2 2 1 1 2 2 1 2 2 2
## 463 464 465 466 467 468 469 471 473 474 475 476 477 478 479 480 481 482 483 484
## 2 2 2 2 2 2 1 1 1 2 1 1 2 2 1 2 1 1 2 1
## 485 487 488 489 490 491 492 493 494 495 496 497 498 499 500 502 504 505 506 507
## 1 2 1 1 1 2 2 2 2 2 2 1 2 2 2 2 1 2 1 1
## 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527
## 2 1 1 1 1 2 2 2 2 2 2 1 2 1 1 2 2 1 2 1
## 528 529 531 532 533 534 535 536 537 538 539 540 541 542 544 545 546 547 549 550
## 1 2 2 1 2 1 1 1 2 1 1 1 1 2 1 2 1 1 1 1
## 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
## 1 1 1 2 1 1 2 1 2 2 1 1 2 2 1 1 1 1 1 1
## 571 573 574 575 576 577 578 580 581 583 584 585 586 587 588 589 590 591 592 593
## 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1
## 594 595 596 597 603 605 611 613 614 617 619 621 624 625 629 632 638 640 644 646
## 1 2 1 1 2 2 1 2 2 1 1 1 1 2 1 1 1 1 2 2
## 651 652 653 655 657 658 659 663 668 672 674 678 679 684 685 687 688 689 693 699
## 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 701 702 705 710 711 712 714 716 717 719 721 722 723 724 725 726 727 728 730 731
## 1 1 1 1 1 1 1 1 2 1 1 1 2 2 1 1 1 1 2 2
## 732 733 735 736 737 738 739 740 741 742 743 745 746 748 749 752 753 755 756 758
## 1 1 2 1 1 2 1 1 1 1 1 2 1 2 1 1 1 1 1 2
## 760 761 762 764 766 767 768 771 773 774 775 776 780 784 785 787 788 792 794 795
## 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 2 1
## 796 797 799 800 801 802 803 804 805 806 807 809 811 812 813 814 816 819 820 821
## 1 2 1 2 1 2 1 1 2 2 2 1 1 1 2 1 1 1 2 1
## 822 824 825 826 827 829 837 838 839 842 843 844 847 848 849 851 852 853 855 856
## 1 2 1 2 1 2 1 2 2 1 1 1 1 1 1 2 1 2 1 1
## 860 861 864 867 869 871 872 873 874 875 877 883 887 890 891 892 893 894 896 897
## 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
## 898 899 900 901 903 904 905 907 908 909 910 911 912 913 914 915 916 918 920
## 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1Para visualizar los clústers podemos usar la función clusplot. Vemos la agrupación de 2 a 8 clúster
Observamos estos gráficos podemos comprobar que para k=5 u k =6 se observan grupos diferenciados.
También podemos visualizar el resultado del proceso de agregación con el siguiente código para el caso de 2 clusters
plot(totalData[y_cluster2==1,],col='blue', xlim=c(min(totalData[,1]), max(totalData[,1])), ylim=c(min(totalData[,2]), max(totalData[,2])))
points(totalData[y_cluster2==2,],col='red')
para 5:
plot(totalData[y_cluster5==1,],col='blue', xlim=c(min(totalData[,1]), max(totalData[,1])), ylim=c(min(totalData[,2]), max(totalData[,2])))
points(totalData[y_cluster5==2,],col='red')
points(totalData[y_cluster5==3,],col='green')
points(totalData[y_cluster5==4,],col='yellow')
points(totalData[y_cluster5==5,],col='black')
y para 6:
plot(totalData[y_cluster6==1,],col='blue', xlim=c(min(totalData[,1]), max(totalData[,1])), ylim=c(min(totalData[,2]), max(totalData[,2])))
points(totalData[y_cluster6==2,],col='red')
points(totalData[y_cluster6==3,],col='green')
points(totalData[y_cluster6==4,],col='black')
points(totalData[y_cluster6==5,],col='yellow')
points(totalData[y_cluster6==6,],col='purple')
Ahora vamos a evaluar la calidad del proceso de agregación. Para ello usaremos la función silhouette que calcula la silueta de cada muestra
d <- daisy(totalData)
sk2 <- silhouette(y_cluster2, d)
sk3 <- silhouette(y_cluster3, d)
sk4 <- silhouette(y_cluster4, d)
sk5 <- silhouette(y_cluster5, d)
sk6 <- silhouette(y_cluster6, d)
sk7 <- silhouette(y_cluster7, d)
sk8 <- silhouette(y_cluster8, d)La función silhouette devuelve para cada muestra, el cluster dónde ha sido asignado, el cluster vecino y el valor de la silueta. Por lo tanto, calculando la media de la tercera columna podemos obtener una estimación de la calidad del agrupamiento
mean(sk2[,3])## [1] 0.3702318mean(sk3[,3])## [1] 0.3763959mean(sk4[,3])## [1] 0.3537122mean(sk5[,3])## [1] 0.3342239mean(sk6[,3])## [1] 0.3625661mean(sk7[,3])## [1] 0.3451166mean(sk8[,3])## [1] 0.3422829Como se puede comprobar, k = 2 seguida de k = 6 es la mejor opción.
Como inicialmente no conocemos el número óptimo de clústers, probamos con varios valores
d <- daisy(totalData)
resultados <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
fit <- kmeans(totalData, i)
y_cluster <- fit$cluster
sk <- silhouette(y_cluster, d)
resultados[i] <- mean(sk[,3])
}Mostramos en un gráfica los valores de las siluetas media de cada prueba para comprobar que número de clústers es el mejor
plot(2:10,resultados[2:10],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="Silueta")
Verificamos que el mejor valor que se obtiene es k=2, k=3.
Otro forma de evaluar cual es el mejor número de clústers es considerar el mejor modelo, aquel que ofrece la menor suma de los cuadrados de las distancias de los puntos de cada grupo con respecto a su centro (withinss), con la mayor separación entre centros de grupos (betweenss). Como se puede comprobar es una idea conceptualmente similar a la silueta. Una manera común de hacer la selección del número de clústers consiste en aplicar el método elbow (codo), que no es más que la selección del número de clústers en base a la inspección de la gráfica que se obtiene al iterar con el mismo conjunto de datos para distintos valores del número de clústers. Se seleccionará el valor que se encuentra en el “codo” de la curva
resultados <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
fit <- kmeans(totalData, i)
resultados[i] <- fit$tot.withinss
}
plot(2:10,resultados[2:10],type="o",col="blue",pch=0,xlab="Número de clusters",ylab="tot.tot.withinss")
En este caso el número óptimo de clústers es 5 o 6, ya que es cuando la curva comienza a estabilizarse.
También se puede usar la función kmeansruns del paquete fpc que ejecuta el algoritmo kmeans con un conjunto de valores, para después seleccionar el valor del número de clústers que mejor funcione de acuerdo a dos criterios: la silueta media (“asw”) y Calinski-Harabasz (“ch”).
library(fpc)
fit_ch <- kmeansruns(totalData, krange = 1:10, criterion = "ch")
fit_asw <- kmeansruns(totalData, krange = 1:10, criterion = "asw") Podemos comprobar el valor con el que se ha obtenido el mejor resultado y también mostrar el resultado obtenido para todos los valores de k usando ambos criterios
fit_ch$bestk## [1] 1fit_asw$bestk## [1] 1plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio Calinski-Harabasz")
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio silueta media")
Los resultados son muy parecidos a los que hemos obtenido anteriormente. Con el criterio de la silueta media se obtienen 2,3 y 6 clústers y con el Calinski-Harabasz se obtienen 3 y 6.
Como se ha comprobado, conocer el número óptimo de clústers no es un problema fácil. Tampoco lo es la evaluación de los modelos de agregación.
Tras las pruebas que hemos realizado para obtener el número óptimo de clusters, hemos encontrado que el número de clusters varía según el método entre 3,2,5, o 6. Vamos a estudiar los resultados encontrados entre 2 y 6 clusters.
A continuación mostramos visualmente los clusters encontrados suponiendo que hay 2, 3 y 6 clusters. Hay que tener en cuenta que nos es posible mostrar los clusters en un espacio de 4 dimensiones. Por lo tanto, mostramos los clusters entre pares de características, ademas en este modelo solo usamos dos variables:
cl3 <- kmeans(totalData, 2)
with(totalData, pairs(totalData, col=c(1:2)[cl3$cluster])) 
cl3 <- kmeans(totalData, 3)
with(totalData, pairs(totalData, col=c(1:3)[cl3$cluster])) 
cl3 <- kmeans(totalData, 4)
with(totalData, pairs(totalData, col=c(1:4)[cl3$cluster])) 
cl3 <- kmeans(totalData, 5)
with(totalData, pairs(totalData, col=c(1:5)[cl3$cluster])) 
cl3 <- kmeans(totalData, 6)
with(totalData, pairs(totalData, col=c(1:6)[cl3$cluster])) 
if (!require('factoextra')) install.packages('factoextra'); library('factoextra')library(ggpubr)
library(factoextra)# Compute k-means with k = 5
set.seed(123)
res.km <- kmeans(scale(totalData[, -5]), 6, nstart = 25)
# K-means clusters showing the group of each individuals
res.km$cluster## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 4 3 5 5 5 5 5 2 2 2 2 2 5 6 6 3 5 2 5
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 3 6 5 5 2 3 5 4 6 1 2 3 2 5 5 5 1 4 2 2
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 4 5 6 6 5 3 3 2 2 2 3 3 3 5 2 1 5 3 3 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 2 2 2 5 1 2 2 6 6 6 6 5 1 3 5 6 3 2 5 4
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 3 2 5 6 5 5 2 1 5 2 3 6 2 5 5 5 3 6 5 5
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 5 5 5 3 3 3 5 2 3 3 2 3 5 2 1 2 5 5 2 2
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 2 6 3 1 5 5 6 3 5 5 3 3 5 5 5 5 2 1 3 5
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 5 6 5 2 3 3 4 5 5 3 6 3 3 6 1 1 5 5 5 3
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 3 5 3 3 5 5 5 5 5 3 4 4 6 2 2 4 3 1 5 5
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 5 2 6 6 6 5 5 4 6 2 5 2 2 1 1 3 4 2 5 4
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
## 3 6 6 2 3 2 2 2 5 6 5 5 5 6 3 5 3 2 1 5
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## 5 3 3 1 3 5 3 6 1 3 5 4 3 1 6 1 1 3 5 5
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
## 3 5 5 2 1 1 3 1 5 2 3 4 1 3 5 5 3 2 6 3
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
## 3 2 6 5 2 2 5 2 6 2 2 6 2 3 5 6 2 2 6 2
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
## 3 5 2 5 6 3 6 3 5 5 6 5 3 2 3 5 4 2 3 2
## 301 302 303 304 305 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
## 1 5 5 5 5 6 3 3 5 5 5 5 6 3 5 2 5 6 5 5
## 322 323 324 325 326 327 328 329 330 332 333 334 336 337 339 340 341 342 343 344
## 5 6 3 2 5 1 2 5 5 2 2 5 5 3 4 5 6 3 2 2
## 345 346 347 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
## 5 2 5 3 5 2 5 5 3 3 3 3 3 2 3 6 5 3 5 5
## 366 367 368 370 371 372 373 374 375 377 378 380 381 382 383 384 386 387 388 389
## 5 1 6 2 1 3 2 6 5 2 2 3 2 1 5 3 3 3 4 2
## 391 392 393 396 397 398 399 400 402 403 404 407 408 409 410 413 414 415 416 417
## 3 6 2 1 1 5 1 2 5 2 3 5 2 3 5 3 1 5 5 4
## 418 419 420 421 422 423 424 425 426 427 429 430 431 432 433 434 436 437 439 440
## 2 2 6 1 5 4 3 2 1 3 2 2 1 2 6 5 2 3 5 3
## 441 442 443 444 445 446 447 449 450 451 452 453 454 455 456 457 459 460 461 462
## 1 2 2 3 3 3 5 2 4 3 3 3 1 4 6 3 2 3 3 2
## 463 464 465 466 467 468 469 471 473 474 475 476 477 478 479 480 481 482 483 484
## 3 2 2 3 3 3 1 1 2 3 1 4 2 2 1 5 1 2 2 1
## 485 487 488 489 490 491 492 493 494 495 496 497 498 499 500 502 504 505 506 507
## 4 3 1 1 4 2 3 3 5 3 5 2 3 5 3 3 2 3 1 1
## 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527
## 2 4 1 1 4 6 5 5 2 2 6 1 2 1 2 2 2 2 3 4
## 528 529 531 532 533 534 535 536 537 538 539 540 541 542 544 545 546 547 549 550
## 1 3 3 2 3 2 1 1 3 4 4 4 1 3 1 3 1 2 6 2
## 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
## 4 4 4 2 2 1 6 1 5 5 1 1 2 2 1 1 2 4 2 4
## 571 573 574 575 576 577 578 580 581 583 584 585 586 587 588 589 590 591 592 593
## 1 6 4 3 4 4 1 2 4 1 2 2 1 6 1 4 2 2 2 4
## 594 595 596 597 603 605 611 613 614 617 619 621 624 625 629 632 638 640 644 646
## 1 6 4 1 3 3 4 3 2 1 4 4 1 3 1 1 1 1 2 3
## 651 652 653 655 657 658 659 663 668 672 674 678 679 684 685 687 688 689 693 699
## 1 1 1 4 1 3 1 2 1 1 1 1 4 1 1 4 4 6 4 4
## 701 702 705 710 711 712 714 716 717 719 721 722 723 724 725 726 727 728 730 731
## 4 6 1 4 2 2 2 1 2 4 1 2 2 2 1 1 4 6 6 2
## 732 733 735 736 737 738 739 740 741 742 743 745 746 748 749 752 753 755 756 758
## 2 1 2 1 4 3 1 4 1 1 1 5 2 2 1 1 4 1 2 3
## 760 761 762 764 766 767 768 771 773 774 775 776 780 784 785 787 788 792 794 795
## 4 1 5 1 1 1 1 1 4 5 1 2 4 4 1 2 1 2 2 1
## 796 797 799 800 801 802 803 804 805 806 807 809 811 812 813 814 816 819 820 821
## 4 2 4 2 1 2 4 2 2 2 2 4 1 1 3 1 4 1 5 2
## 822 824 825 826 827 829 837 838 839 842 843 844 847 848 849 851 852 853 855 856
## 4 3 2 3 2 6 1 3 2 2 1 1 1 4 2 5 1 3 4 1
## 860 861 864 867 869 871 872 873 874 875 877 883 887 890 891 892 893 894 896 897
## 4 4 1 2 2 2 1 1 4 1 2 4 1 4 1 5 2 1 1 6
## 898 899 900 901 903 904 905 907 908 909 910 911 912 913 914 915 916 918 920
## 2 4 2 1 3 1 1 1 4 4 2 1 4 2 6 2 5 1 1Para finalizar realizaremos una representación gráfica de los 5 grupos, en función de la puntuación de gasto y los ingresos.
fviz_cluster(res.km, data = totalData[, -5],
palette = c("#2E9FDF", "purple", "#E7B800","red","green","black"),
geom = "point",
ellipse.type = "convex",
ggtheme = theme_bw())
A continuación vamos a utilizar dos métodos distintos del kmeans: el método Hierarchical clustering y el método K-Medois clustering(PAM).
El método Hierarchical clustering es un método de análisis de grupos puntuales, el cual busca construir una jerarquía de grupos.
Dentro de este método encontramos dos subtipos:
Agglomerative clustering (bottom-up): Empieza considerando que cada objeto forma un grupo por sí mismo y a partir de ahí , evalúan las distancias entre grupos y crean los diferentes grupos finales por aglomeración.
Y Divisive clustering (top-down):Este es un acercamiento descendente: todas las observaciones comienzan en un grupo, y se realizan divisiones mientras uno baja en la jerarquía.
En este caso utilizaremos el método de bottom-up.
En primer lugar,escalaremos los datos ya que este método esta basados en medidas de la distancia
Data_escalados <- scale(totalData)Utilizaremos los datos previamente escalados para calcular la matriz de distancias euclidianas:
mat_dist <- dist(Data_escalados , method = "euclidean")En tercer lugar calculamos los dendrogramas con linkage complete y average
de_complete <- hclust(d = mat_dist, method = "complete")
de_verage <- hclust(d = mat_dist, method = "average")cor(x = mat_dist, cophenetic(de_complete))## [1] 0.5641461Con la variable trestbps este valor era de 0.564, vemos que si eliminamos esta variable aumenta considerablemente.
cor(x = mat_dist, cophenetic(de_verage))## [1] 0.6684025Con la variable edad este valor era de 0.6683198, vemos que si eliminamos esta variable aumenta considerablemente.
Las observaciones serán más similares conforme nos acercamos a 1.
Calcularemos la distancia cophenetic del HC con el fin de conocer cual método representa mejor la similitud entre las observaciones, el método complete o average.
hc_euclidea_completo <- hclust(d = dist(Data_escalados, method = "euclidean"),
method = "complete")hc_euclidea_media <- hclust(d = dist(Data_escalados, method = "euclidean"),
method = "average")distancia <- dist(x = Data_escalados, method = "euclidean")distanciac<- cor(x = distancia, cophenetic(hc_euclidea_completo))
distanciaa<- cor(x = distancia, cophenetic(hc_euclidea_media))
distanciac## [1] 0.5641461distanciaa## [1] 0.6684025LA distancia average es mayor por lo que lo ideal sería elegir este método.
Vamos a representar gráficamente el modelo para k = 2,3 y 6, y para linkage: complete y average.
k = 2
fviz_dend(x = hc_euclidea_completo, k = 2, cex = 0.6) +
geom_hline(yintercept = 5.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage complete, K=2")
fviz_dend(x = hc_euclidea_media, k = 2, cex = 0.6) +
geom_hline(yintercept = 5.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage average, K=2")
k = 3
fviz_dend(x = hc_euclidea_completo, k = 3, cex = 0.6) +
geom_hline(yintercept = 3.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage complete, K=3")
fviz_dend(x = hc_euclidea_media, k = 3, cex = 0.6) +
geom_hline(yintercept = 3.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage average, K=3")
k = 6
fviz_dend(x = hc_euclidea_completo, k = 6, cex = 0.6) +
geom_hline(yintercept = 3.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage complete, K=6")
fviz_dend(x = hc_euclidea_media, k = 6, cex = 0.6) +
geom_hline(yintercept = 3.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Linkage average, K=6")
Observamos para finalizar el gráfico Hierarchical clustering + Proyección PCA, para k = 2 y 6.
fviz_cluster(object = list(data=Data_escalados, cluster=cutree(hc_euclidea_completo, k=2)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídea, Linkage complete, K=2") +
theme_bw() +
theme(legend.position = "bottom")## Warning: ggrepel: 623 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_cluster(object = list(data=Data_escalados, cluster=cutree(hc_euclidea_media, k=2)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídea, Linkage average, K=2") +
theme_bw() +
theme(legend.position = "bottom")## Warning: ggrepel: 623 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_cluster(object = list(data=Data_escalados, cluster=cutree(hc_euclidea_completo, k=6)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídea, Linkage complete, K=6") +
theme_bw() +
theme(legend.position = "bottom")## Warning: ggrepel: 629 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_cluster(object = list(data=Data_escalados, cluster=cutree(hc_euclidea_media, k=6)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídea, Linkage average, K=6") +
theme_bw() +
theme(legend.position = "bottom")## Warning: ggrepel: 629 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
El segundo método que vamos a utilizar es K-medoids, este algoritmo se relaciona con el algoritmo k-means y medoifshift. K-medoid es una técnica clásica de particionado de grupos que divide los datos conformados por n objetos en k grupos (con k conocido de antemano).
Es más robusto ante el ruido y a partes aisladas que k-means porque minimiza una suma de disimilaridades (entre pares de puntos) en vez de una suma de distancias euclidianas cuadradas.
Debido a estas características es utilizado en los modelos en los que existen outliers, por ello en este método es útil usar como medida de similitud la distancia Manhatan.
Uno de los algoritmos más utilizados en k-memoids es PAM (Partitioning Around Medoids).
Utilizaremos los datos escalados del modelo anterior y realizamos una repesantación gráfica utilizando la matriz de distancias manhattan
fviz_nbclust(x = Data_escalados, FUNcluster = pam, method = "wss", k.max = 15,
diss = dist(Data_escalados, method = "manhattan"))
Podemos decir que la curva se estabiliza entre 5 y 9 clusters.
Probaremos con 6 clusters ya que erá la mejor opción de los modelos anteriores.
set.seed(123)
pam_clusters <- pam(x = Data_escalados, k = 6, metric = "manhattan")
pam_clusters## Medoids:
## ID Trestbps_norm Thalach_norm
## 480 456 0.3939011 0.82289849
## 423 406 0.9513368 -0.72409206
## 454 433 -0.7209703 -1.11083970
## 27 27 -0.7209703 1.28699566
## 900 723 -0.1635346 0.04940322
## 531 503 -1.2784060 0.43615085
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 2 3 4 4 4 1 4 5 1 1 1 5 4 1 1 6 1 5 4
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 6 1 4 4 5 4 4 2 1 3 1 4 1 1 4 1 3 2 5 2
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 2 1 1 1 4 6 6 2 1 5 6 5 6 4 5 3 1 6 5 5
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 5 1 5 1 3 5 1 1 2 2 1 4 3 6 4 1 5 5 4 2
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 6 5 1 1 4 1 1 3 1 5 6 2 5 4 1 1 6 1 1 4
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 4 4 1 6 5 6 1 5 5 5 1 5 4 5 3 5 1 1 5 5
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 5 1 6 2 1 4 2 6 4 4 6 6 4 1 4 1 2 3 3 4
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 1 1 4 5 6 6 2 4 4 6 1 3 6 2 3 3 1 4 1 6
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 6 4 6 3 4 1 1 1 5 6 2 2 2 1 2 2 6 3 1 4
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 4 5 1 2 1 1 4 2 1 1 1 2 5 3 3 6 2 1 4 2
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
## 6 1 1 5 6 1 5 2 5 1 4 4 4 1 6 4 4 1 5 1
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## 4 6 4 3 4 4 6 1 3 6 1 2 3 3 1 3 3 5 1 4
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
## 6 4 1 5 3 3 6 3 4 5 6 2 3 4 4 4 6 2 2 5
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
## 6 1 1 4 2 5 5 5 1 5 5 2 2 3 1 1 1 1 1 5
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
## 6 4 5 4 1 6 2 5 1 4 1 1 5 1 5 4 2 2 6 5
## 301 302 303 304 305 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
## 3 4 1 4 4 1 6 6 4 5 4 4 1 6 4 1 4 1 4 4
## 322 323 324 325 326 327 328 329 330 332 333 334 336 337 339 340 341 342 343 344
## 4 1 3 5 4 3 5 4 1 1 2 4 4 5 2 4 1 6 5 5
## 345 346 347 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
## 1 5 4 6 4 5 4 4 5 6 6 6 4 2 6 2 1 6 4 4
## 366 367 368 370 371 372 373 374 375 377 378 380 381 382 383 384 386 387 388 389
## 4 3 1 5 3 5 5 1 4 5 2 5 1 3 1 6 6 6 2 5
## 391 392 393 396 397 398 399 400 402 403 404 407 408 409 410 413 414 415 416 417
## 6 1 2 3 3 4 3 2 1 2 6 4 5 4 1 6 3 4 1 2
## 418 419 420 421 422 423 424 425 426 427 429 430 431 432 433 434 436 437 439 440
## 5 2 1 3 1 2 5 5 3 6 5 5 3 5 1 1 5 3 1 5
## 441 442 443 444 445 446 447 449 450 451 452 453 454 455 456 457 459 460 461 462
## 3 5 1 5 5 6 5 2 2 5 6 4 3 2 1 5 2 6 4 5
## 463 464 465 466 467 468 469 471 473 474 475 476 477 478 479 480 481 482 483 484
## 5 1 1 6 5 5 2 3 5 5 3 2 5 5 3 1 3 5 1 5
## 485 487 488 489 490 491 492 493 494 495 496 497 498 499 500 502 504 505 506 507
## 2 5 5 3 2 1 6 6 1 6 4 5 6 4 6 4 2 6 3 3
## 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527
## 1 2 3 3 2 1 4 4 5 1 1 5 5 3 5 1 5 5 6 2
## 528 529 531 532 533 534 535 536 537 538 539 540 541 542 544 545 546 547 549 550
## 2 3 6 2 4 5 3 3 6 2 2 2 5 6 3 6 3 5 2 2
## 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
## 2 2 2 5 5 3 1 3 1 1 3 5 5 1 3 3 2 2 2 2
## 571 573 574 575 576 577 578 580 581 583 584 585 586 587 588 589 590 591 592 593
## 2 2 2 6 2 2 3 2 2 3 5 5 3 1 3 2 5 2 5 2
## 594 595 596 597 603 605 611 613 614 617 619 621 624 625 629 632 638 640 644 646
## 3 1 2 3 6 3 2 6 5 3 2 2 3 4 3 2 3 3 5 4
## 651 652 653 655 657 658 659 663 668 672 674 678 679 684 685 687 688 689 693 699
## 3 2 3 2 3 6 3 2 3 3 3 3 2 3 3 2 2 2 2 2
## 701 702 705 710 711 712 714 716 717 719 721 722 723 724 725 726 727 728 730 731
## 2 2 3 2 2 2 5 3 1 2 2 5 5 1 3 3 2 2 1 5
## 732 733 735 736 737 738 739 740 741 742 743 745 746 748 749 752 753 755 756 758
## 2 3 5 3 2 6 3 2 3 3 3 4 5 5 3 3 2 3 2 5
## 760 761 762 764 766 767 768 771 773 774 775 776 780 784 785 787 788 792 794 795
## 2 3 4 3 3 3 3 3 2 4 3 5 2 2 3 2 3 1 5 3
## 796 797 799 800 801 802 803 804 805 806 807 809 811 812 813 814 816 819 820 821
## 2 5 2 5 3 1 2 5 5 5 5 2 3 3 6 3 2 3 1 2
## 822 824 825 826 827 829 837 838 839 842 843 844 847 848 849 851 852 853 855 856
## 2 3 5 5 2 1 3 6 5 2 3 3 3 2 2 1 2 6 2 3
## 860 861 864 867 869 871 872 873 874 875 877 883 887 890 891 892 893 894 896 897
## 2 2 3 1 1 2 3 3 2 3 2 2 3 2 3 4 2 2 3 1
## 898 899 900 901 903 904 905 907 908 909 910 911 912 913 914 915 916 918 920
## 5 2 5 3 3 3 3 3 2 2 5 3 2 2 2 5 5 3 3
## Objective function:
## build swap
## 0.7379849 0.7241316
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"Observamos gráficamente los resultados del método PAM
fviz_cluster(object = pam_clusters, data = Data_escalados, ellipse.type = "t",
repel = TRUE) +
theme_bw() +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")## Warning: ggrepel: 655 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Para finalizar vamos a comparar los distintos métodos utilizados
Observando los distintos modelos utilizados podemos dar por buenos los resultados del k = 6, encontramos 6 grupos bastante diferenciados de pacientes con unas características peculiares.
Estudiaremos las características de los grupos diferenciados para comprender el modelo resultante:
Para estudiar los grupos que se han creado tendremos en cuenta lo siguiente:
Dividiremos el Trestbps en 3 partes:
De 0.0 a 0.2 Trestbps normalizado: lo llamaremos Trestbps_norm_bajo
De 0.2 a 0.6 Trestbps normalizado: lo llamaremos Trestbps_norm_medio
De 0.6 a 1 Trestbps normalizado: lo llamaremos Trestbps_norm_alto
Haremos lo mismo para Thalach_norm, en este caso en :
De 0.0 a 0.4 Thalach normalizado: lo llamaremos Thalach_norm_bajo
De 0.4 a 0.7 Thalach normalizado: lo llamaremos Thalachs_norm_medio
De 0.7 a 1 Thalach normalizado: lo llamaremos Trestbps_norm_alto
De esta manera agruparemos de la siguiente manera. 6 tipos de pacientes diferenciados por su *presión arterial en reposo y frecuencia máxima alcanzada**:
Pacientes tipo 1 : Trestbps_norm_bajo + Thalach_norm_bajo
Pacientes tipo 2 : Trestbps_norm_bajo y Thalachs_norm_medio
Pacientes tipo 3 : Trestbps_norm_bajo y Thalach_norm_bajo
Pacientes tipo 4 : Trestbps_norm_medio y Thalach_norm_medio
Pacientes tipo 5 : Trestbps_norm_alto y Thalach_norm_alto
Pacientes tipo 6 : Trestbps_norm_alto y Thalach_norm_bajo
Estos grupos servirán de utilidad y se podrían utilizar para futuros análisis o investigaciones médicas.
Previamente en la Prac anterior, modificamos el modelo para poder aplicar el estudio de reglas de asociación
Para poder llevar a cabo el algoritmo Apriori deberemos tener una estructura de los datos de tipo transacción.
Cargamos la librería necesaria:
library(arules)Cargamos y transformamos el dataframe para la asociación a una transacción
diagnostico_paciente <- read.transactions(file = "data/pacientes_diagnosticos.csv",
format = "single",
sep = ",",
header = TRUE,
cols = c("ID_Paciente", "Value"),
rm.duplicates = TRUE)Observamos el resumen de la transacción:
inspect(head(diagnostico_paciente,3))## items transactionID
## [1] {angina_típica,
## pendiente_ST_descendente,
## resultados_electro_hipertrofia,
## sin_enfermedad_cardiaca,
## talasemia_defecto_fijo,
## vasos_coloreados_0} 1
## [2] {asintomatico,
## con_enfermedad_cardiaca,
## pendiente_ST_descendente,
## resultados_electro_hipertrofia,
## talasemia_defecto_reversible,
## vasos_coloreados_0} 10
## [3] {asintomatico,
## pendiente_ST_ascendente,
## resultados_electro_hipertrofia,
## sin_enfermedad_cardiaca,
## talasemia_normal,
## vasos_coloreados_0} 100summary(diagnostico_paciente)## transactions as itemMatrix in sparse format with
## 739 rows (elements/itemsets/transactions) and
## 20 columns (items) and a density of 0.2631935
##
## most frequent items:
## resultados_electro_normal desconocido asintomatico
## 444 440 392
## con_enfermedad_cardiaca sin_enfermedad_cardiaca (Other)
## 382 357 1875
##
## element (itemset/transaction) length distribution:
## sizes
## 4 5 6
## 178 188 373
##
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.000 5.000 6.000 5.264 6.000 6.000
##
## includes extended item information - examples:
## labels
## 1 angina_atípica
## 2 angina_típica
## 3 asintomatico
##
## includes extended transaction information - examples:
## transactionID
## 1 1
## 2 10
## 3 100De la tabla anterior obtenemos la siguiente información:
- Disponemos de 739 usuarios
- Disponemos de 6 diagnósticos.
- El diagnóstico más común es resultados_electro_normal (444), seguido por desconocido (440).
- Cada paciente tiene 6 tipos de diagnósticos, hay 739 pacientes con seis diagnósticos cada uno.En el gráfico inferior observamos los diagnósticos más frecuentes:
itemFrequencyPlot(diagnostico_paciente,topN=30,type="absolute")
Los diagnósticos más frecuentes corresponden a la situación de normalidad de cada estudio
En este momento ya podemos utilizar el algoritmo apriori para estudiar las reglas de asociación:
El algoritmo apriori esta compuesto por tres elementos:
Support:hace referencia a la popularidad de un determinado diagnóstico.
Confidence:probabilidad de que también un paciente tenga un diagnóstico ‘x’ si tiene otro ‘y’.
Lift: En este caso hace referencia al aumento de tener un diagnóstico ‘x’ cuando tenemos un diagnóstico ‘y’.
reglas_diagnosticos <- apriori(diagnostico_paciente, parameter = list(support = 0.01, confidence = 0.5))## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## 0.5 0.1 1 none FALSE TRUE 5 0.01 1
## maxlen target ext
## 10 rules TRUE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 7
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[20 item(s), 739 transaction(s)] done [0.00s].
## sorting and recoding items ... [20 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 6 done [0.00s].
## writing ... [1293 rule(s)] done [0.00s].
## creating S4 object ... done [0.00s].Se han creado 1923 reglas. Observaremos las 5 primeras ordenadas por elementos
Vamos a proceder a ordenarlas según su Support:
support <- inspect(head(sort(reglas_diagnosticos, by = "support"), 5))## lhs rhs support
## [1] {} => {resultados_electro_normal} 0.6008119
## [2] {} => {desconocido} 0.5953992
## [3] {} => {asintomatico} 0.5304465
## [4] {} => {con_enfermedad_cardiaca} 0.5169147
## [5] {con_enfermedad_cardiaca} => {asintomatico} 0.4032476
## confidence coverage lift count
## [1] 0.6008119 1.0000000 1.000000 444
## [2] 0.5953992 1.0000000 1.000000 440
## [3] 0.5304465 1.0000000 1.000000 392
## [4] 0.5169147 1.0000000 1.000000 382
## [5] 0.7801047 0.5169147 1.470657 298Confirmamos que la popularidad del diagnostico hacen referencias a diagnósticos positivos, aquellos donde no se han encontrado problemas en los pacientes y diagnósticos desconocidos.
confidence <- inspect(head(sort(reglas_diagnosticos, by = "confidence"), 5))## lhs rhs support confidence coverage lift count
## [1] {angina_típica,
## talasemia_defecto_reversible} => {vasos_coloreados_0} 0.01082544 1 0.01082544 4.105556 8
## [2] {talasemia_defecto_reversible,
## vasos_coloreados_2} => {con_enfermedad_cardiaca} 0.02706360 1 0.02706360 1.934555 20
## [3] {angina_atípica,
## resultados_electro_anomalia} => {desconocido} 0.02165088 1 0.02165088 1.679545 16
## [4] {pendiente_ST_ascendente,
## resultados_electro_anomalia} => {desconocido} 0.01759134 1 0.01759134 1.679545 13
## [5] {resultados_electro_hipertrofia,
## talasemia_defecto_reversible,
## vasos_coloreados_3} => {con_enfermedad_cardiaca} 0.01082544 1 0.01082544 1.934555 8Observamos las reglas por Confidence donde destacamos la regla 2:
SI el paciente sufre talasemia_defecto_reversible y vasos_coloreados_2 será probable que sufra una enfermedad cardíaca. Y la regla 5: resultados_electro_hipertrofia, talasemia_defecto_reversible,vasos_coloreados_3, llevaran a sufrir una enfermedad cardíaca.
Y por Lift:
lift <- inspect(head(sort(reglas_diagnosticos, by = "lift"), 5))## lhs rhs support confidence coverage lift count
## [1] {asintomatico,
## con_enfermedad_cardiaca,
## resultados_electro_hipertrofia,
## talasemia_normal} => {vasos_coloreados_1} 0.01082544 0.5 0.02165088 5.684615 8
## [2] {asintomatico,
## pendiente_ST_plano,
## resultados_electro_normal,
## vasos_coloreados_1} => {talasemia_defecto_reversible} 0.01217862 1.0 0.01217862 4.247126 9
## [3] {asintomatico,
## con_enfermedad_cardiaca,
## pendiente_ST_plano,
## resultados_electro_normal,
## vasos_coloreados_1} => {talasemia_defecto_reversible} 0.01082544 1.0 0.01082544 4.247126 8
## [4] {angina_típica,
## talasemia_defecto_reversible} => {vasos_coloreados_0} 0.01082544 1.0 0.01082544 4.105556 8
## [5] {resultados_electro_hipertrofia,
## sin_enfermedad_cardiaca,
## talasemia_defecto_reversible} => {vasos_coloreados_0} 0.01217862 1.0 0.01217862 4.105556 9Observaremos estas reglas de forma visual:
support_5 <- head(sort(reglas_diagnosticos, by = "support"), 5)
confidence_5 <- head(sort(reglas_diagnosticos, by = "confidence"), 5)
lift_5 <- head(sort(reglas_diagnosticos, by = "lift"), 5)library(arulesViz)plot(support_5, method = "graph", engine = "htmlwidget" )plot(confidence_5, method = "graph", engine = "htmlwidget")plot(lift_5, method = "graph", engine = "htmlwidget")Podemos comprobar que las reglas de asociación tienen sentido ya que para el método support, los diagnósticos más habituales como ya hemos señalado son los diagósticos de normalidad.
Para el método cofidence y lift deberíamos valorar con un experto en cardiología el sentido de estas reglas, ya que al ser un modelo creado por nosotros podría llevar a alguna regla de asociación extraña.
Observamos el dataframe preparado anteriormente para los algoritmos de clasificación:
head (dfc)## Age Gender CP Trestbps RestECG Exang Oldpeak Slope Diagnosis FCMax_Thalac
## 1 63 1 1 145 2 0 2.3 3 0 1
## 2 67 1 4 160 2 1 1.5 2 1 1
## 3 67 1 4 120 2 1 2.6 2 1 1
## 4 37 1 3 130 0 0 3.5 3 0 0
## 5 41 0 2 130 2 0 1.4 1 0 1
## 6 56 1 2 120 0 0 0.8 1 0 0
## edad_discretizada chol_discretizado trestbps_discretizado
## 1 2 1 3
## 2 3 2 3
## 3 3 1 0
## 4 0 2 1
## 5 1 1 1
## 6 2 1 0
## oldpeak_discretizado
## 1 2
## 2 1
## 3 2
## 4 2
## 5 1
## 6 1str(dfc)## 'data.frame': 739 obs. of 14 variables:
## $ Age : num 63 67 67 37 41 56 62 57 63 53 ...
## $ Gender : num 1 1 1 1 0 1 0 0 1 1 ...
## $ CP : num 1 4 4 3 2 2 4 4 4 4 ...
## $ Trestbps : num 145 160 120 130 130 120 140 120 130 140 ...
## $ RestECG : num 2 2 2 0 2 0 2 0 2 2 ...
## $ Exang : num 0 1 1 0 0 0 0 1 0 1 ...
## $ Oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ Slope : num 3 2 2 3 1 1 3 1 2 3 ...
## $ Diagnosis : num 0 1 1 0 0 0 1 0 1 1 ...
## $ FCMax_Thalac : num 1 1 1 0 1 0 0 0 1 1 ...
## $ edad_discretizada : num 2 3 3 0 1 2 2 2 2 1 ...
## $ chol_discretizado : num 1 2 1 2 1 1 2 2 2 1 ...
## $ trestbps_discretizado: num 3 3 0 1 1 0 2 0 1 2 ...
## $ oldpeak_discretizado : num 2 1 2 2 1 1 2 1 1 2 ...vamos a ordenar las columnas por orden de importancia tomando como referencia el análisis PCA de la primera parte del modelo
dfc <- dfc[,c(14,7,1,11,6,8,13,3,4,5,2,10,12,9)] Mostramos el dataframe ordenado que utilizaremos para llevar a cabo los arboles de decisión
head(dfc)## oldpeak_discretizado Oldpeak Age edad_discretizada Exang Slope
## 1 2 2.3 63 2 0 3
## 2 1 1.5 67 3 1 2
## 3 2 2.6 67 3 1 2
## 4 2 3.5 37 0 0 3
## 5 1 1.4 41 1 0 1
## 6 1 0.8 56 2 0 1
## trestbps_discretizado CP Trestbps RestECG Gender FCMax_Thalac
## 1 3 1 145 2 1 1
## 2 3 4 160 2 1 1
## 3 0 4 120 2 1 1
## 4 1 3 130 0 1 0
## 5 1 2 130 2 0 1
## 6 0 2 120 0 1 0
## chol_discretizado Diagnosis
## 1 1 0
## 2 2 1
## 3 1 1
## 4 2 0
## 5 1 0
## 6 1 0Nos interesa desordenar los datos para poder elegir registros aleatorios. Guardaremos los datos con el nuevo nombre como dr1 para el árbol de decisión con todas las variables anteriores.
set.seed(1)
dr1 <- dfc[sample(nrow(dfc)),]Para la futura evaluación del árbol de decisión, es necesario dividir el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba. El conjunto de entrenamiento es el subconjunto del conjunto original de datos utilizado para construir un primer modelo; y el conjunto de prueba, el subconjunto del conjunto original de datos utilizado para evaluar la calidad del modelo.
Lo más correcto será utilizar un conjunto de datos diferente del que utilizamos para construir el árbol, es decir, un conjunto diferente del de entrenamiento. No hay ninguna proporción fijada con respecto al número relativo de componentes de cada subconjunto, pero la más utilizada acostumbra a ser 2/3 para el conjunto de entrenamiento y 1/3, para el conjunto de prueba.
La variable por la que clasificaremos es Diagnosis, para clasificar si existe riesgo de enfermedad cardíaca o no, y que se sitúa en nuestro dataframe en la última columna, la 14. De esta forma, tendremos un conjunto de datos para el entrenamiento y uno para la validación
head(dr1)## oldpeak_discretizado Oldpeak Age edad_discretizada Exang Slope
## 820 2 3.0 63 2 0 2
## 129 0 0.0 44 1 0 1
## 537 1 2.0 48 1 1 3
## 496 2 3.0 36 0 0 2
## 299 1 1.2 45 1 0 2
## 270 0 0.0 42 1 0 1
## trestbps_discretizado CP Trestbps RestECG Gender FCMax_Thalac
## 820 1 3 130 1 1 0
## 129 0 2 120 0 1 1
## 537 1 4 122 1 1 1
## 496 0 2 120 0 1 1
## 299 0 1 110 0 1 1
## 270 1 3 130 0 1 1
## chol_discretizado Diagnosis
## 820 0 0
## 129 1 0
## 537 2 1
## 496 2 1
## 299 2 1
## 270 0 0set.seed(666)
y <- dr1[,14]
X <- dr1[,1:13] Dividimos el dataset en dos partes diferenciadas, train y test. Podemos crear directamente un rango utilizando el parámetro split_prop
split_prop <- 3
indexes = sample(1:nrow(dr1), size=floor(((split_prop-1)/split_prop)*nrow(dr1)))
trainx<-X[indexes,]
trainy<-y[indexes]
testx<-X[-indexes,]
testy<-y[-indexes]Después de una extracción aleatoria de casos es altamente recomendable efectuar un análisis de datos mínimo para asegurarnos de no obtener clasificadores sesgados por los valores que contiene cada muestra. En este caso, verificaremos que la proporción del defaults es más o menos constante en los dos conjuntos:
summary(trainx);## oldpeak_discretizado Oldpeak Age edad_discretizada
## Min. :0.0000 Min. :-1.0000 Min. :28.00 Min. :0.000
## 1st Qu.:0.0000 1st Qu.: 0.0000 1st Qu.:46.00 1st Qu.:1.000
## Median :1.0000 Median : 0.5000 Median :54.00 Median :1.000
## Mean :0.6768 Mean : 0.8957 Mean :53.12 Mean :1.419
## 3rd Qu.:1.0000 3rd Qu.: 1.5000 3rd Qu.:59.25 3rd Qu.:2.000
## Max. :3.0000 Max. : 5.6000 Max. :77.00 Max. :3.000
## Exang Slope trestbps_discretizado CP
## Min. :0.0000 Min. : 1.000 Min. :0.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.: 2.000 1st Qu.:0.000 1st Qu.:3.000
## Median :0.0000 Median : 2.000 Median :1.000 Median :4.000
## Mean :0.4045 Mean : 4.089 Mean :1.384 Mean :3.268
## 3rd Qu.:1.0000 3rd Qu.:10.000 3rd Qu.:2.000 3rd Qu.:4.000
## Max. :1.0000 Max. :10.000 Max. :4.000 Max. :4.000
## Trestbps RestECG Gender FCMax_Thalac
## Min. : 92 Min. :0.0000 Min. :0.0000 Min. :0.000
## 1st Qu.:120 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.000
## Median :130 Median :0.0000 Median :1.0000 Median :1.000
## Mean :133 Mean :0.6057 Mean :0.7561 Mean :0.876
## 3rd Qu.:140 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :200 Max. :2.0000 Max. :1.0000 Max. :1.000
## chol_discretizado
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :1.154
## 3rd Qu.:2.000
## Max. :2.000summary(trainy)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.5142 1.0000 1.0000summary(testx)## oldpeak_discretizado Oldpeak Age edad_discretizada
## Min. :0.0000 Min. :0.0000 Min. :29.00 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:47.00 1st Qu.:1.000
## Median :1.0000 Median :0.4000 Median :54.00 Median :1.000
## Mean :0.6883 Mean :0.8891 Mean :53.05 Mean :1.389
## 3rd Qu.:1.0000 3rd Qu.:1.5000 3rd Qu.:60.00 3rd Qu.:2.000
## Max. :4.0000 Max. :6.2000 Max. :77.00 Max. :3.000
## Exang Slope trestbps_discretizado CP
## Min. :0.0000 Min. : 1.000 Min. :0.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.: 1.500 1st Qu.:0.000 1st Qu.:2.000
## Median :0.0000 Median : 2.000 Median :1.000 Median :3.000
## Mean :0.3927 Mean : 4.101 Mean :1.356 Mean :3.146
## 3rd Qu.:1.0000 3rd Qu.:10.000 3rd Qu.:2.500 3rd Qu.:4.000
## Max. :1.0000 Max. :10.000 Max. :4.000 Max. :4.000
## Trestbps RestECG Gender FCMax_Thalac
## Min. : 96.0 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:120.0 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.0000
## Median :130.0 Median :0.0000 Median :1.0000 Median :1.0000
## Mean :132.9 Mean :0.6964 Mean :0.7814 Mean :0.8543
## 3rd Qu.:141.0 3rd Qu.:2.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :190.0 Max. :2.0000 Max. :1.0000 Max. :1.0000
## chol_discretizado
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :1.194
## 3rd Qu.:2.000
## Max. :2.000summary(testy)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.5223 1.0000 1.0000Verificamos que no hay diferencias graves entre train y test que puedan sesgar las conclusiones.
Se crea el árbol de decisión usando los datos de entrenamiento:
if(!require(ggplot2)){
install.packages('ggplot2', repos='http://cran.us.r-project.org')
library(ggplot2)
}
if(!require(ggpubr)){
install.packages('ggpubr', repos='http://cran.us.r-project.org')
library(ggpubr)
}
if(!require(grid)){
install.packages('grid', repos='http://cran.us.r-project.org')
library(grid)
}## Loading required package: gridif(!require(gridExtra)){
install.packages('gridExtra', repos='http://cran.us.r-project.org')
library(gridExtra)
}## Loading required package: gridExtra##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combineif(!require(C50)){
install.packages('C50', repos='http://cran.us.r-project.org')
library(C50)
}## Loading required package: C50En un primer análisis vamos a realizar el árbol sin poda, noGlobalPruning = FALSE y lo compararemos con el árbol podado en el apartado siguiente
trainy = as.factor(trainy)
model <- C50::C5.0(trainx, trainy,rules=TRUE,control = C5.0Control(noGlobalPruning = FALSE) )
summary(model)##
## Call:
## C5.0.default(x = trainx, y = trainy, rules = TRUE, control
## = C5.0Control(noGlobalPruning = FALSE))
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 15 23:27:31 2021
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 492 cases (14 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (14, lift 1.9)
## oldpeak_discretizado <= 1
## edad_discretizada <= 2
## Slope <= 1
## CP > 2
## CP <= 3
## chol_discretizado <= 1
## -> class 0 [0.938]
##
## Rule 2: (14, lift 1.9)
## Oldpeak <= 1.6
## edad_discretizada > 0
## Exang <= 0
## trestbps_discretizado <= 0
## Trestbps > 110
## RestECG <= 0
## chol_discretizado > 1
## -> class 0 [0.938]
##
## Rule 3: (29/1, lift 1.9)
## Slope > 1
## CP <= 3
## Gender <= 0
## chol_discretizado <= 1
## -> class 0 [0.935]
##
## Rule 4: (43/2, lift 1.9)
## Exang <= 0
## Slope <= 2
## trestbps_discretizado <= 2
## Gender <= 0
## -> class 0 [0.933]
##
## Rule 5: (109/8, lift 1.9)
## Exang <= 0
## trestbps_discretizado <= 1
## CP <= 3
## -> class 0 [0.919]
##
## Rule 6: (10, lift 1.9)
## oldpeak_discretizado > 0
## oldpeak_discretizado <= 1
## Exang <= 0
## Slope <= 1
## trestbps_discretizado > 1
## chol_discretizado > 1
## -> class 0 [0.917]
##
## Rule 7: (80/6, lift 1.9)
## edad_discretizada <= 2
## Exang <= 0
## Slope > 2
## CP <= 3
## -> class 0 [0.915]
##
## Rule 8: (14/1, lift 1.8)
## Exang <= 0
## trestbps_discretizado > 0
## trestbps_discretizado <= 1
## RestECG <= 0
## chol_discretizado > 1
## -> class 0 [0.875]
##
## Rule 9: (29/3, lift 1.8)
## CP <= 3
## Trestbps <= 115
## -> class 0 [0.871]
##
## Rule 10: (5, lift 1.8)
## Exang > 0
## Slope <= 1
## CP <= 3
## -> class 0 [0.857]
##
## Rule 11: (5, lift 1.8)
## Oldpeak <= 0.9
## Exang > 0
## CP > 3
## Gender <= 0
## -> class 0 [0.857]
##
## Rule 12: (34/6, lift 1.7)
## Exang <= 0
## Slope > 3
## chol_discretizado <= 0
## -> class 0 [0.806]
##
## Rule 13: (184/44, lift 1.6)
## Oldpeak <= 1.6
## Exang <= 0
## chol_discretizado > 0
## -> class 0 [0.758]
##
## Rule 14: (23, lift 1.9)
## oldpeak_discretizado > 1
## Slope <= 2
## trestbps_discretizado > 1
## -> class 1 [0.960]
##
## Rule 15: (16, lift 1.8)
## Oldpeak > 0.4
## Oldpeak <= 1.6
## CP > 3
## RestECG > 1
## Gender > 0
## -> class 1 [0.944]
##
## Rule 16: (57/3, lift 1.8)
## Slope <= 3
## CP > 3
## Gender > 0
## chol_discretizado <= 0
## -> class 1 [0.932]
##
## Rule 17: (72/5, lift 1.8)
## Oldpeak > 1.6
## CP > 3
## Gender > 0
## -> class 1 [0.919]
##
## Rule 18: (22/2, lift 1.7)
## edad_discretizada > 2
## trestbps_discretizado > 1
## Gender > 0
## -> class 1 [0.875]
##
## Rule 19: (37/5, lift 1.6)
## edad_discretizada <= 2
## Slope <= 2
## trestbps_discretizado > 1
## Gender > 0
## chol_discretizado > 1
## -> class 1 [0.846]
##
## Rule 20: (4, lift 1.6)
## Exang <= 0
## trestbps_discretizado > 2
## CP > 3
## Gender <= 0
## -> class 1 [0.833]
##
## Rule 21: (68/11, lift 1.6)
## trestbps_discretizado <= 1
## CP > 3
## RestECG <= 0
## Gender > 0
## -> class 1 [0.829]
##
## Rule 22: (199/38, lift 1.6)
## Exang > 0
## -> class 1 [0.806]
##
## Rule 23: (89/22, lift 1.5)
## Age > 55
## trestbps_discretizado > 1
## Gender > 0
## -> class 1 [0.747]
##
## Default class: 1
##
##
## Evaluation on training data (492 cases):
##
## Rules
## ----------------
## No Errors
##
## 23 55(11.2%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 208 31 (a): class 0
## 24 229 (b): class 1
##
##
## Attribute usage:
##
## 91.46% Exang
## 67.28% CP
## 66.06% trestbps_discretizado
## 65.04% chol_discretizado
## 60.77% Gender
## 55.08% Oldpeak
## 54.47% Slope
## 33.13% edad_discretizada
## 21.75% RestECG
## 18.09% Age
## 9.55% oldpeak_discretizado
## 8.74% Trestbps
##
##
## Time: 0.0 secsErrors muestra el número y porcentaje de casos mal clasificados en el subconjunto de entrenamiento. El árbol obtenido clasifica erróneamente 58 de los 492 casos dados (2/3) deL total de observaciones del dataframe, correspondiente a una tasa de error de 11.2% .
Más abajo podemos observar los atributos usados de más a menos porcentaje de uso. Destacando en el top 5:
91.46% Exang
67.28% CP
66.06% trestbps_discretizado
65.04% chol_discretizado
60.77% Gender
55.08% Oldpeak
A partir del árbol de decisión que hemos modelado, vamos a describir las reglas que superan la probabilidad de acierto del 90%:
Vamos a observar aquellas reglas con un nivel de precisión mayor al 90% sobre la posibilidad de sufrir una enfermedad cardíaca:
Regla 1: Si oldpeak_discretizado > 1, Slope <= 2, trestbps_discretizado > 1 hay un 96% de posibilidades sufrir una enfermedad cardíaca.
Regla 2: Si Oldpeak > 0.4, Oldpeak <= 1.6, CP > 3, RestECG > 1, Gender > 0class 1 hay un 94.4% de posibilidades sufrir una enfermedad cardíaca.
Regla 3: Si Slope <= 3, CP > 3, Gender > 0 (hombre), chol_discretizado <= 0 hay un 93.2% de posibilidades sufrir una enfermedad cardíaca.
Regla 4: Si Oldpeak > 1.6, CP > 3, Gender > 0 (hombre), hay un 91.9% de posibilidades sufrir una enfermedad cardíaca.
A continuación, mostramos el árbol obtenido.
model <- C50::C5.0(trainx, trainy)
plot(model)
Una vez tenemos el modelo, podemos comprobar su calidad prediciendo la clase para los datos de prueba que nos hemos reservado al principio.
predicted_model <- predict( model, testx, type="class" )
print(sprintf("La precisión del modelo es: %.4f %%",100*sum(predicted_model == testy) / length(predicted_model)))## [1] "La precisión del modelo es: 78.5425 %"Cuando hay pocas clases, la calidad de la predicción se puede analizar mediante una matriz de confusión que identifica los tipos de errores cometidos.
mat_conf<-table(testy,Predicted=predicted_model)
mat_conf## Predicted
## testy 0 1
## 0 98 20
## 1 33 96Además, tenemos a nuestra disposición el paquete gmodels para obtener información más completa:
if(!require(gmodels)){
install.packages('gmodels', repos='http://cran.us.r-project.org')
library(gmodels)
}## Loading required package: gmodelsCrossTable(testy, predicted_model,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 247
##
##
## | Prediction
## Reality | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 98 | 20 | 118 |
## | 0.397 | 0.081 | |
## -------------|-----------|-----------|-----------|
## 1 | 33 | 96 | 129 |
## | 0.134 | 0.389 | |
## -------------|-----------|-----------|-----------|
## Column Total | 131 | 116 | 247 |
## -------------|-----------|-----------|-----------|
##
## Fijándonos en la tabla superior, hemos clasificado 53 registros incorrectamente de 247, el 78.5425 % por ciento.
Otra manera de calcular el porcentaje de registros correctamente clasificados es usando la matriz de confusión:
porcentaje_correct<-100 * sum(diag(mat_conf)) / sum(mat_conf)
print(sprintf("El %% de registros correctamente clasificados es: %.4f %%",porcentaje_correct))## [1] "El % de registros correctamente clasificados es: 78.5425 %"Hemos obtenido una tasa de acierto del 78.5425 %, vamos a intentar mejorarla en los siguientes modelos.
Vamos a crear un modelo de árbol con un número de ramas menor, intentando alcanzar una visualización más eficiente. El algoritmo C50::C5.0 realiza la poda automáticamente si añadimos en el parámetro control, noGlobalPruning = TRUE.
trainy = as.factor(trainy)
modelo3 <- C50::C5.0(trainx, trainy,rules=TRUE,control = C5.0Control(noGlobalPruning = TRUE) )predicted_model3 <- predict( modelo3, testx, type="class" )
print(sprintf("La precisión del modelo3 es: %.4f %%",100*sum(predicted_model3 == testy) / length(predicted_model3)))## [1] "La precisión del modelo3 es: 76.9231 %"La crosstable del modelo 3 es:
CrossTable(testy, predicted_model3,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 247
##
##
## | Prediction
## Reality | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 95 | 23 | 118 |
## | 0.385 | 0.093 | |
## -------------|-----------|-----------|-----------|
## 1 | 34 | 95 | 129 |
## | 0.138 | 0.385 | |
## -------------|-----------|-----------|-----------|
## Column Total | 129 | 118 | 247 |
## -------------|-----------|-----------|-----------|
##
## La precisión se ha reducido por la poda pasando de 78.5425 % a 76.9231 %".
Utilizaremos el enfoque de boosting basado en el trabajo de Rob Schapire and Yoav Freund (1999) y gracias al parámetro trials, del algoritmo C50::C5.0
modelo2 <- C50::C5.0(trainx, trainy, trials = 100)predicted_model2 <- predict( modelo2, testx, type="class" )
print(sprintf("La precisión del modelo 2 es: %.4f %%",100*sum(predicted_model2 == testy) / length(predicted_model2)))## [1] "La precisión del modelo 2 es: 78.9474 %"Con el número máximo de intentos que nos permite el modelo, 100 intentos, observamos como se aumenta la precisión del modelo, llegando al 78.9474 %.
La crosstable del modelo 2 será;
CrossTable(testy, predicted_model2,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 247
##
##
## | Prediction
## Reality | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 99 | 19 | 118 |
## | 0.401 | 0.077 | |
## -------------|-----------|-----------|-----------|
## 1 | 33 | 96 | 129 |
## | 0.134 | 0.389 | |
## -------------|-----------|-----------|-----------|
## Column Total | 132 | 115 | 247 |
## -------------|-----------|-----------|-----------|
##
## En este momento vamos a utilizar otro método para crear nuestro árbol de decisión.
Utilizaremos dos bibliotecas diferentes, una de ellas es rpart, para crear el árbol de decisión y rpart.plot para representar visualmente nuestro árbol de decisión.
library("rpart")
library("rpart.plot") modelo4 <- rpart(formula = trainy ~ ., data = trainx, method = "class" )En el primer árbol Observamos las tasas de clasificación para cada nodo, expresada como el número de clasificaciones correctas y el número de observaciones en el nodo. En el segundo gráfico observamos la probabilidad por cada clase.
rpart.plot(modelo4, type=2,extra = 2, under = TRUE, faclen=5,cex=.55) 
rpart.plot(modelo4, type=2,extra = 8, under = TRUE, faclen=5,cex=.55) 
En el primer árbol observamos las tasas de clasificación para cada nodo, expresada como el número de clasificaciones correctas y el número de observaciones en el nodo.
En el segundo gráfico observamos la probabilidad por cada clase. Realizamos la predicción con este nuevo método
Para algunas clases, clasificamos correctamente mas del 80% de las veces.
predicted_model4 <- predict( modelo4, testx, type="class" )
print(sprintf("La precisión del modelo4 es: %.4f %%",100*sum(predicted_model4 == testy) / length(predicted_model4)))## [1] "La precisión del modelo4 es: 74.8988 %"La crosstable del modelo 4 será;
CrossTable(testy, predicted_model4,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 247
##
##
## | Prediction
## Reality | 0 | 1 | Row Total |
## -------------|-----------|-----------|-----------|
## 0 | 95 | 23 | 118 |
## | 0.385 | 0.093 | |
## -------------|-----------|-----------|-----------|
## 1 | 39 | 90 | 129 |
## | 0.158 | 0.364 | |
## -------------|-----------|-----------|-----------|
## Column Total | 134 | 113 | 247 |
## -------------|-----------|-----------|-----------|
##
## set.seed(666)
y1 <- dr1[,14]
X1 <- dr1[,1:13] Dividimos el dataset en dos partes diferenciadas, train y test. Podemos crear directamente un rango utilizando el parámetro split_prop
split_prop <- 3
indexes = sample(1:nrow(dr1), size=floor(((split_prop-1)/split_prop)*nrow(dr1)))
trainx1<-X1[indexes,]
trainy1<-y1[indexes]
testx1<-X1[-indexes,]
testy1<-y1[-indexes]trainy1 = as.factor(trainy1)
modelo5 <- C50::C5.0(trainx1, trainy1,rules=TRUE, trials = 100 )predicted_model5 <- predict( modelo5, testx1, type="class" )
print(sprintf("La precisión del modelo5 es: %.4f %%",100*sum(predicted_model5 == testy1) / length(predicted_model5)))## [1] "La precisión del modelo5 es: 80.1619 %"Reduciendo las variables no nos aumenta la precisión del modelo máxima que hemos alcanzado
Hemos conseguido obtener como precisión máxima un 80 % , además de las reglas comentadas anteriormente que nos ofrecen características clínicas de aquellos pacientes con más probabilidad de sufrir una enfermedad cardíaca.
Los mejores resultados obtenidos son los siguientes: 80.1619 %, que corresponderá al modelo con boosting, con reglas y sin poda. Seguido por el modelo de boosting con poda 78.947 % y el modelo C50 normal sin poda 78.5425 %.
D considerar que este modelo estaba pre-diseñado para la realización de un modelo de clasificación, donde a través de varios dataframes, con datos de pacientes relacionados con las enfermedades cardíacas, clasificar estos pacientes entre aquellos que tendrán riesgo de sufrir una enfermedad cardíaca y los que no. Por ello, no hay ninguna limitación aparente para los modelos de clasificación.
En referencia a los modelos de clustering, al tratarse de datos sobre pacientes, podríamos decir que dichas características ya estarán agrupadas por si mismas, por ello no tendría mucha utilidad realizar este análisis para este tipo de datos.
Para el modelo de reglas de asociación, si que sería interesante realizar comparaciones entre diagnósticos, ya que nos ofrecería información importante de sobre el comportamiento de los pacientes, pero necesitaríamos muchos más datos para que el modelo fuera consistente.
Por otro lado, para que el modelo fuera más eficiente, deberíamos añadir más variables y parámetros al modelo. El número de causas externas que afectan a los pacientes en distintas enfermedades es elevado, a veces causas como la climatología, el periodo temporal, etc, pueden ser muy influyentes en este tipo de enfermedades.
Sociedad Española de Cardiología (SEC):
Enfermedades cardiovasculares:
INE (Defunciones por causa de muerte):
Para la PRAC_2 he utilizado como fuentes mis PECs anteriores de la asignatura